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ABSTRACT 


The computer code G3DXL computes time domain scattered fields and 
surface currents and charges induced by a driving function on and within a 
complex scattering object which may be perfectly conducting or a lossy 
dielectric. This is accomplished by modeling the object with cells within a 
three-dimensional, rectangular problem space, enforcing the appropriate 
boundary conditions and differencing Maxwell's equations in time. In the 
present version of the program, the driving function can be either the field 
radiated by a lightning strike or a direct lightning strike. The scattering 
object is the F-106 B aircraft. 

The program has been successfully executed on the CDC 203 at NASA 
Langley where it is currently operational. 
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SECTION 1 


SCOPE AND OBJECTIVES 


This document presents information concerning the programming, 
maintenance and operation of the computer code G3DXL. This information is 
intended to facilitate 1) effective use of the program as is and 2) conversion 
and modification of the code for other test items and excitation sources. 

The theoretical basis for the code is discussed, the numerical 
implementation of the theory is presented, and the limitations are enumerated. 
This discussion is followed by a presentation of the program architecture, 
including flow charts, global variable definitions, input and output 
requirements. A discussion of each individual subroutine is given which 
includes a subroutine description, local variable definition (where needed), 
flow chart, and a listing of the subroutine. The reader is provided in this 
manual with an understanding of the theoretical basis of the code and its 
practical implementation. 



SECTION 2 


PROBLEM DEFINITION 


2.1 THE PROBLEM 

G3DXL addresses the problem of calculating in the time domain the 
currents and charges induced on and within a composite aircraft, along with the 
resulting scattering fields, by an excitation source, in this case a direct 
lightning strike or the radiated lightning field. The scatterer may be 
imbedded in a homogeneous media (an aircraft in flight, for example) or may be 
located over a lossy ground plane (an aircraft over a runway, for example). 

2.2 THE SOLUTION 

G3DXL finds the induced currents and charges on the scatterer along 
with the scattered fields by differencing Maxwell's equations in time using an 
algorithm first described by Yee [1]. This algorithm is applied to a problem 
space with a mesh of dimensions of 29x29x29. Even larger problem spaces on the 
order of 50x50x50 are possible on certain large computers, such as the DEC VAX 
11/780. The problem of echoes from the finite problem space's outer boundaries 
is circumvented by the application of Merewether's radiating boundary condition 
[2], which is based on the fact that, far from the scatterer, the scattered 
fields must behave as f(0,<J>)g(t-r/c)/r. Here, f(6,<J>) represents an arbitrary 
variation in 6 and <t>, g(t-r/c) expresses the retarded time behavior of the 
fields and the 1/r describes the far field dependence on r. It is the 1/r 
dependence that is of important here, as it allows fields near the outer 
boundary to be extrapolated to the outer boundary, thereby affecting the outer 
radiating boundary condition. 
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2.3 Numerical Implementation 


The three-dimensional, finite-difference formalism employed by G3DXL 
separates the fields into incident and scattered fields and solves directly for 
the scattered fields, using Maxwell's equations. The incident field is used 
only to establish the value of the scattered field throughout the scatterer's 
volume, where the fields are related by what will be referred to here as a 
volume boundary condition. 

The formalism (Appendix B gives more details) starts with Maxwell's 
equations for the total fields: 


» * * - -4 

-► 3E -► 

V x H = e~ + oE. 

Subtracting the Maxwell equations for the incident field propagating in free 
space, i.e. the free space form of the Maxwell equations: 


V x E 1 = -uoff- 


V x H 1 = e 0 |f- 


yields the generalized equations: 


-+53 -H 

3H n +s , . 3H 

u— = -v x e - (y-Mo)^- 


9E 3 +S n \ ^ 

+oE = VxH - oE - (e-e 0 )^— 
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The generalized equations, as they stand, are completely general and 
may be applied over all space. In free space, the terms with (y-y Q ), (e-e Q ) 
and o as multipliers vanish yielding the free space form of the equations. On 
a perfectly conducting surface, where o=°° the generalized equations reduce to: 



where T refers to the tangential components of the E-field on the perfectly 
conducting surface. 


Prior work with finite difference codes typically used the free field 
formulation with the perfectly conducting boundary condition. G3DXL uses the 
generalized equations throughout the volume of a lossy dielectric scatterer, 
the free field equations in free space and the perfectly conducting boundary 
condition for a perfect conductor. This approach was deemed more efficient 
than the use of the generalized equations throughout all space. 


The actual conversion of the differential equations into a finite 
difference form will be presented here for a perfectly conducting scatterer 
about which the free space Maxwell equations for the scattered field hold, 
i.e. , 


V x E 3 


9H S 

‘at 


„ 3 s 3E S jjs 
V x » = ^ + aE 

but are written in a finite-difference form using two point central 
differencing in a rectangular coordinate system. For example: 


3H scat 3H scat 

__z 1 

3y 3z 


3E 


scat 


3t 


+ oE 


scat 

x 
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becomes 


( H^ cat ) " ( I ,J,K) - (H^ Cat ) ( I , J-1 ,K) 

Y(J) - Y( J-1 ) 

(H aCat ) (I ,J ,K) - (H aCat ) ( I , J , K-1 ) 

Z ( K) - Z(K-1) 

„ n+1/2 „ n-1/2 

(E aCat ) (I , J ,K) - (E aCat ) ( I , J , K) 


, n+1/2 . n-1/2 

(E a ° at ) (I ,J ,K) + (E aCat ) (I,J,K) 


where At, Y ( J ) - Y(J~1), etc. are discrete time and spatial increments used in 
the calculation. The superscript n is the time index and I,J,K label the 
nodes in the three-dimensional mesh of the problem space. Similar expressions 
can be derived for the remaining five equations, involving the remaining 
components. Note that the H-field components at times (n+1)At can be evaluated 
from their own value at the earlier time nAt and from the E-field components at 
the earlier times (n+1/2)At. Similarly the E-field components at times 
(n+1/2)At are evaluated from the H-field components at times nAt and from their 
own values at earlier times (n-1/2)At. This results in a fully explicit and 
automatically central differencing scheme that reduces the truncation or round- 
off error [3.^1 . 

The difference equations are applied, at successive discrete time 
steps, to a cell space; in this case, a three-dimensional, rectangular cell 
space. That is, the equations are assumed to hold for all the cells of the 
space. A single cell and the field components associated with it are shown in 
Figure 1. The size of the individual cell determines the time step At. For 
stability, At must satisfy the Courant stability condition, which may be 
expressed as [5]. 
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Figure 1. Location of the Six Field Evaluation 
Points in a Typical Cell 




At = min . „( (X (1+1) - X (I)) 2 (eu) 1 
I , J ,K o o 

+ (Y (J+1) - Y (J)) _2 (ey) _1 + (Z (K+1) - Z (K) )" 2 ( ey ) _1 ) _1 /2 
o o o o 

If £ or p are inhomogenous , their dependence on (I,J,K) must be included in 
evaluating the minimal time steps of the above equation. If the time step is 
greater than At, the predicted responses grow without bound. Note that At is 
approximately the transit time across the cell. The cell size also determines 
the frequency range of valid data. When a quarter wavelength of the frequency 
component of interest is smaller than the cell dimensions, the answer returned 
is in error. 

Two boundary conditions must be applied to this cell space. The 
first defines the object for which the scattered fields are to be calculated. 
We shall only consider a perfect conductor here for simplicity. For each cell 
face that defines a portion of aircraft surface, the total tangential E-field 
vanishes so that 


scat 

E tan 


(t) 


inc 

tan 


(t) 


At t=0 , E* nc (t) is assumed zero and so E SCa ^(t) is also zero. At later times 

inc 

the source or driving field E (t) appears at the scatterer's surface with 

tan 

scat 

non-zero value determining E (t) and, through the difference equations, 

tan 

determining the scattered fields throughout the cell space. 


The second boundary condition is on the outer surface of the cell 

space where the radiation condition is imposed [2]. This boundary condition 

scat 

requires the scattered E-field to behave as E‘ (t) = f(6,$) g(t-r/c)/r at the 
outer surface, the same field behavior as occurs far from the scatter. By 
imposing this behavior on the surface of the cell space, a finite cell space 
can be made to approximate an infinite cell space. In effect, the outermost E- 
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fields are made to behave nearly the same as they would have had the cell space 
continued on beyond the outermost cells indefinitely, leaving the inner cells 
seemingly imbedded in an infinite cell space. As a rule of thumb, this 
approximation is successful when the cell space is at least twice the scatter's 
dimensions in each direction. 


As an example of this procedure, let J be the maximum index in the 

m n+1 

y-coordinate for the cell space. The radiation condition applied to E in 

z 

the y-direction is: 


R <I,J-1,K) 

E" (I,J ,K) = 
z m 


2R d,j ,k) * 


nr 


+ 2(1-6 2 )E n (I,J -1,K) + 6 (0 +1)E n+1 (I,J -1 ,K) ] 

zy z m zy zy z m 


where 

R (I,J ,K) = \/x (I ) 2 + Y (J ) 2 + Z (K ) 2 
z m Vo o m o 

R (I, J ,K) - R (I, J “1 , K) 

. z m z m 

zy cAt 

and X (I), Y (J), Z (K) are coordinate values for the points considered, c is 
o o o 

the velocity of light. As can be seen, the fields at the outer boundary are 
found from interpolating, in time, the fields that are one cell in from the 
outer boundary. The times are the present program cycle, and the two prior 
program cycles. This interpolation is based on a constant cell size. If an 
expanding mesh is used, a somewhat more complex expression must be employed 
[5]. 

Similar conditions (there are eleven more) can be written down for 

E , E „ and E at the six surfaces of the box. 
x’ y’ z 
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By assigning appropriate constitutive parameters (e, p, and o) to the 
first few tiers of cells at the base of the cell space, a ground can be modeled 
from which the scattered fields are partially reflected. Since the radiation 
condition is applied to all six faces, only the desired ground reflection will 
be present; the fields transmitted into the ground will not be reflected back 
from the outer surface of the cell space. 

The cell size used need not be uniform. By enlarging the cell size 
outside the immediate region of the aircraft, the outer boundary can be further 
removed from the aircraft. (Expansion factors much greater than 1.3 are not 
recommended. The reason is that severe discontinuity effects can be seen.) 
Reflections from the outer boundary, that occur because the radiation condition 
is imposed at moderate distance R from the origin, are thereby lessened. 

The surface current density and surface charge density at any 
position on the aircraft are found from 


J 


n x 


,-scat — inc ^ , 

(H + H ) and p 


E scat 

normal 


+ 


E inc 

normal 


respectively, where n is the outward normal of the scatterer's surface. The 
values of the scattered fields at the position of interest are found from the 
adjacent nodal point values using a combination of interpolation and 
extrapolation. These time-domain responses can be transformed to yield 
frequency domain data over a broad frequency range. 


4 
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Program Archi tectural Description 


2.4 


G3DXL is designed to predict exterior and interior responses on 
complex scattering objects, such as aircraft. To achieve the needed degree of 
spatial resolution for interior coupling predictions the expansion technique is 
employed. A preliminary run with a coarsely modeled aircraft (Figure 1-a) is 
used to obtain the response about a portion of the aircraft. This response 
imposed on the outer boundary of the problem space allows a second run to be 
made in which only the portion of the aircraft is modeled (Figure 2), but in 
finer detail and with the entire aircraft response preserved. The use of the 
expansion technique requires what is, in effect, two separate runs. 

The scattering object can also be a perfect conductor, a restriction 
that some prior finite difference codes, such as THREDE [5,6], had to adhere 
to; or it could be a lossy dielectric, as in a composite aircraft. This last 
and more general case is treated using the generalized finite difference 
technique as embodied in G3DXL. A thin plate algorithm (Appendix A) even 
allows composite panels to be treated. The generalized finite difference 
technique and thin plate algorithm do not require more runs, but each run using 
them is more complex. Such a run must be able to treat perfectly conducting 
solids and surfaces or plates and must be able to treat lossy dielectric solids 
and plates. 



Unexpanded Model 
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Bui khead 




2.5 Program Execution 


The following steps are a concise guide to the execution of the code. 

1) Always maintain two original backup copies: the first, a direct 

access file resident in the computer’s memory and the second, a card deck, the 
same as the file, resident with the user in charge of maintaining the code. 

2) From the original file in the computer, make a new file under a 
new name. It is this file that will be modified (using the XEDIT capability) 
to model the correct problem. 

3) When the new file is modified and running properly, it is prudent 
to make a copy of the file. This copy serves as a secondary backup. If 
further modifications are made and they result in a nonrunning program, one can 
always return to the secondary backup and try again. 

4) Modifications to the original code, as embodied in the original 
backup, can be made only after the selection of an excitation source and an 
interaction object. 

5) The excitation source can be either a radiated lightning field or 
a direct lightning strike. To select one or the other, the parameter IDLS is 
set equal to 0 for a radiated lightning field or 1 for a direct lightning 
strike. IDLS appears as a data statement in subroutine SETUP. 

6) The excitation source specifications are set as follows. For a 
radiated lightning field a double exponential waveform A(e at -e is assumed. 
The parameters A, a and 8 are represented by AMP, ALPHA and BETA in the code. 
They are set in subroutine SETUP. Current values are AMP = 5.92E4, ALPHA * 
H.08E6, BETA = 1.0E8. 
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7) The actual calculation of the radiated lightning fields is 

performed in the functions EINCX(I,J,K) HINCZP(XP,YP,ZP) . At present the 

fields are for broadside incidence from above. The electric field points in 
the x direction and the magnetic field points in the z direction. 

8) The interaction object, typically an aircraft, must be modeled 

within the problem space boundaries set by a space of 28 x 28 x 28 cells. 
Further, the outer radiating boundary condition can be expected to work 
reasonably well only if the interaction object occupies, at most, 20 cells 
across the problem space in any direction. Thus, the interaction object can be 

divided by no more than 20 in any linear dimension. If one dimension is L, 

then the cell used in the problem space would have dimension L/20 in this 
direction. A typical aircraft is on the order to 20 meters long. Thus, cells 
approximately 1 meter long are used in the M=1 loop where the entire aircraft 
must be modeled. Generally, cells on the order of lm x 0.5m x 0.5m are used. 
This allows for a reasonably accurate rendition of the aircraft fuselage. 
Expanding the cells where the wings lie is, however, necessary. This allows 
all of the wings to be modeled in the problem space at a sacrifice in detail in 

the region of the wings, where it is often not needed. 

9) For the M=1 loop, when the entire interaction object must be 
modeled, a set of scale drawings or the equivalent are needed. They define the 
interaction object geometry. Additionally, the constitutive parameters (o and 
e) of the material comprising the interaction object are needed along with the 
material thickness. The model is input into the code by: 

o setting the cell size (Ax,Ay,Az) through the parameters DELX, DELY 
and DELZ in subroutine SETUP. 

o setting the expansion rates of the cells (XPANX, XPANY, XPANZ) and 
the indices where cell expansion starts, (IUP, . . .KD0WN) in 
subroutine SETUP. 

o setting the N0PE(I,J,K) array in subroutine BUILD. Nonzero values 
of N0PE(I,J,K) place a material media in the (I,J,K) tn cells. When 
N0PE(I,J,K) = 0 the cell has the properties of free space. The 
option exists to select a perfectly conducting material that fills 
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the cell (NOPE(I,J,K) = 4) or forma a plate on one aide of the cell 
(N0PE(I,J,K) = 1 for the (Y-Z) plane, 2 for the (X-Z) plane and 3 
for the (X-Y) plane). (All aidea are the near aide of the cell aa 
defined along the x,y,z axia that follow the increaae in the 
indicea I, J, and K reapectively. ) The option alao exista to 
replace the perfectly conducting material with a loaay dielectric 
material. The aame approach is uaed aa for the perfectly 
conducting material, except 1 •2,3,4 are replaced by 5,6,7,8 in the 
values that may be aaaigned to N0PE(I,J,K). Celia are aet to 
nonzero values so as to fill the problem with material following 
the geometry of the interaction object and having the same material 
properties (e and a). 

o setting e and o for any loaay dielectric material in the 

interaction model via a data statement for EPS and SIGMA in the 
subroutine SETUP. 

o setting plate thickness TX for lossy dielectric plates in the (Y-Z) 
plane, TY for lossy dielectric plates in (X-Z) plane and TZ for 
lossy dielectric plates in the (X-Y) plane. (Presently EPS, SIGMA, 
TX, TY and TZ do not vary as a function of position). At a large 
cost in memory, but within the capabilities of the Cyber 203, these 
parameters can be input as a function of position, i.e. 

EPS+EPS(I ,J,K) , etc. If this change is required, then an array, 
such as EPS(I,J,K), must be specified in subroutine SETUP and the 
CDC XEDIT capability used to replace EPS by EPS(I,J,K) everywhere, 
including common statements. Related quantities to EPS, such aa 
EPEFFX, EPEFFY, EPEFFZ , EXPON, EXPDX, EXPDY , EXPDZ, XEXP, XEXPX, 
XEXPZ , and to SIGMA, such as RSIGMA , must also be written as three 
dimensional arrays and dimensioned as such or, more appropriately, 
written out explicitly in terms of EPX(I,J,K), SIGMA(I ,J ,K) , 
TX(I,J,K), etc. at the cost of an increased computational burden 
with the savings of a very large amount of memory. 

10) Interaction object models need only be constructed for those 
loops desired, i.e. 

M-1 for unexpanded model where only surface responses are required; 

M«2 for expanded model where only aperture coupling and diffusion 

may be treated; 

11) The desired permutation of the two loops is set by MM. Allowed 
values are MM=1 , 12, where, for example, MM=12 results in the M=1 and 2 loops 
being exercised. MM is set in program DRIVER. 


I 

i 

i 

I 
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12) The number of time steps, N, must also be set in DRIVER. 

Dimension restrictions limit N at present to 2000. Further, the expansion 
factor EXPFAC times N for the M»1 loop should yield N for the M=2 loop. Thus, 
when an expanded run is made, N £ 2000 and N - 2000/EXPFAC (2000/4 = 500 at 
present) for unexpanded runs made in conjunction with an expanded run. 

13) For expanded runs (M*2) the necessary subboundary information on 
the preceding unexpanded run (M=1) is saved via subroutine SAVESB and 
interpolated for use in the expanded run in subroutine 0UTBND. Presently, the 
expansion factor, EXPFAC, is four, so that the space surrounded by the 
subboundary on the unexpanded run is 7 x 7 x 7 cells. The use of interpolation 
only and the positioning of the tangential E field components on the 
subboundary faces requires that an array of 8 x 9 field components be saved on 
each face. As there are two components on each face and six faces, a total of 
864 components are saved in each time step in the array ARAY, which has 
dimensions of (864,500). For different expansion factors the 864 would have to 
change in ARAY, which appears in subroutines SAVESB and 0UTBND. Also, 

CRAY (8, 9) appearing in subroutine 0UTBND for interpolation purposes would have 
to be changed, as would the auxiliary array ARRAY. 

14) The location of the subboundary is set by a data statement in 

INEAR KFAR in subroutines SAVESB and OUTBND . 

15) Interior regions are typically modeled for the M-2 loop where 

SC 3 1 

sufficient resolution is available. Wires are electrically conducting (E‘ 
i nc 

-E ) for M =2. They are manually specified in EBC. 

16) Lightning channels are modeled as exterior wires running from the 
outer surface of the problem space to the scattering object. They are 
electrically conducting for M-1 and 2 and are manually specified in EBC. A 
loop of H-field components drive the electrically conducting lightning channel 
for M-1 and is specified manually in HADV. For M«2 the loop is not presented, 
only the channel. This is because the loop acts as a source and must be 
outside the subboundary volume for the expansion technique to work. 
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17) Test point locations are selected manually in DATASAV. The 
desired locations (up to ZH) are specified by XOBS(NPT), YOBS(NPT), ZOBS(NPT), 
which are the test point locations in meters from the center of the problem 
space. Setting NPLANE(NPT) for each location determines the surface on which 
test point is located or, more precisely, the direction from which the fields 
are extrapolated to the test point. NPLANE=1 selects a test point on a Y-Z 
plane with the fields extrapolated in from left to right; 2 is the same, except 
the extrapolation is from right to left. NPLANE=3 selects a test point on a X- 
Z plane with the fields extrapolated from high to low; is the same except the 
extrapolation is from low to high. NPLANE=5 selects a test point on a X-Y 
plane with the fields extrapolated from back to fore; 6 is the same except the 
extrapolation is from fore to back. A test location on a surface can be for an 
exterior or interior response depending on the selection of NPLANE, which 
determines from what direction the extrapolation is made. THETA(NPT) aligns 
the test point response components, that is for THETA=0, HST0R1 is the axial 
current for the test point on the sides of a cylinder and HST0R2 is the 
circumferential current. When THETA=90, this assignment is reversed. For 
THETA somewhere in between, HST0R1 , for example, is a mix of axial and 
circumferential currents. For most applications THETA can be left at 0. 

18) To obtain current responses for late times economically, the data 
records for the HST0R1 component are extended analytically LMAX * IP time 
steps. IP is the number of time steps taken before response data is saved. 

For example, for IP=2 data is saved every other time step. Both LMAX and IP 
are set via a data statement in subroutine SETUP. For currents on interior 
wires a loop of magnetic field about the wire can be manually specified and 
output as CSTORE. NPLANE must be selected so that each component is stored in 
HST0R1 , which is analytically extended. The, CSTORE is also of extended 
duration and the late time currents are known. 


18 



2.6 


LIMITATIONS 


There are three limitations on the operation of G3DXL which are 
imposed by the mathematical model and the computer facilities. These are 
listed and discussed below. 

1) Run tine is limited to <2000 program cycles for a single 

loop for typical aircraft scattering problems. A program 

cycle is typically on the order of 1 ns resulting in 500 

to 2000 ns of data. Instabilities in the radiation 
boundary condition return a growing oscillation in the 
late-time data that overwhelms the true response. This 
limitation is lessened by maintaining nearly equal cell 
sizes on the outer boundaries (a condition that is 
automatically met when using a constant mesh, but not 
always when using an expanding mesh.) Differences of 
greater than a factor of two should be avoided. 

2) The model is limited in size to be no more than one-half 
the cell space in any direction. If the model is not 
restricted to this size limit, the approximation upon 
which the radiating boundary condition is based, namely 
that the scattered fields behave as far fields at the 
outer boundary, is severely compromised and reflections 
off the outer boundary will be excessively large. 

3) The upper frequency is limited to c/Mi, where J, i3 the 
largest dimension of the cells used to define the 
scatterer. This is equivalent to quarter wavelengths 
equal to the largest cell dimension. At these 
wavelengths, the granularity in the scattering model 
becomes apparent and at shorter wavelengths, i.e., high 
frequencies, invalidates the data. 
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H) Only expansion factors of 2, M, 7, and are allowed, as 
these are the only integral factors possible with a 28- 
on-a-side cell space. Further, the factor of 1 4 is not 
recommended, as it will involve very crude subboundary 
field interpolations. 
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SECTION 3 

PROGRAMMING INFORMATION 


This section provides information directed toward implementation of 
the computer program G3DXL on a particular computer, modification or extension 
to meet particular needs or conversion for a different computer environment. 

3.1 PROBLEM FLOW 

Flow charts showing the overall program structure for G3DXL are given 
in Figures 2a-c. The functions performed by each subroutine shown on the chart 
are described in detail in Section 3.5. The main program inputs a variable 
indicating how frequently to output results; it sets the starting and ending 
time for the computations and also increments time. Most importantly, in 
addition, it calls subprograms SETUP, BUILD, EADV, HADV, SAVESB, and DATASAV. 
The variables input in the main program are described in Table 1 . 

Table 1. Variables Input in the Main Program, DRIVER, of G3DXL 


Variable 

Name 

MM 

IDLS 

TSTART 


Description 

Desired permutation of the four possible configurations 
0=lndirect Radiation 1 “Direct Lightning Strike 
Time (in seconds) at which to start computations 
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Figure 2b. First Pass, Diffusion Calculation 
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Figure 2c. Second Pass, Expanded Diffusion Calculation 
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oooooo o non 


PROGRAM DRIVER( INPUT , OUTPUT-OUTKF ) 

C 

COMMON/ EFIELD/EX( 29 ,29,29) ,EY(29,29,29) ,EZ(29,29,29) 

COMMON/ HFIELD/HX(29 ,29 ,29) ,HY(29 ,29 ,29) ,HZ(29 ,29 ,29) 

COMMON/ GRID/X( 28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX ,NY ,NZ ,HXI ,NY1 ,NZl ,N,M,MQ ,DT,XMU,EPSO ,EPS ,NPTLIM , 

1 NN.NPTS.LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
COMMON/ PERM/ MM 
C 

C LOOP OVER THE DESIRED PERMUTATION OF THE TWO CONFIGURATIONS 
C (UNEXPANDED M=l/ EXPANDED ,M-2) 

C 

C SET MM VALUE TO DETERMINE WHICH PERMUTATION 
C (UNEXPANDED ONL Y, MM- 1/ EXPANDED ,MM=12) 

C 

MM- 12 
C 

PRINT 10 

10 FORMAT(*MM INPUT ERROR*) 

C 

100 CONTINUE 
C 

C UNEXPANDED(M-l) LOOP 

C 

M-l 

C 

C GENERATE PROBLEM SPACE AND INTERACTION OBJECT FOR M-l 
C 

CALL SETUP 
CALL BUILD 

SET TIME LIMITS AND EVERY IP DATA POINT SAVED 

T START-0. 0 
PRINT 111 

111 FORMAT (* BUILD DONE*) 

T-DT/2.+TSTART 
N-0 

CALL EADV 

IF(MM.NE.l) CALL SAVESB 
PRINT 112 

112 FORMAT (* EADV CALLED*) 

DO 130 N-1,48 

ADVANCE TIME 

T-T+DT/2. 

ADVANCE HFIELD 

CALL HADV 
C 
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ooo non oo o non ooo ooo n o o o o ooo ooo 


C ADVANCE TIME 

C 

T=T+DT/2. 

ADVANCE EFIELD 
CALL EADV 

IF(MM.NE.I) CALL SAVESB 
STORE FIELDS 

IF(MOD(N, IP) .EQ.O) CALL DATASAV 


130 CONTINUE 

PRINT 150, T,N 

150 F0RMAT(*1EXIT TIME(M=1)=*E 12. 3* .AFTER CYCLE*I4) 
CALL PRINOUT 
IF(MM.NE.l) GO TO 160 
GO TO 500 

160 CONTINUE 

200 CONTINUE 

EXPANDED (M=2) LOOP 


M=2 

GENERATE PROBLEM SPACE AND INTERACTION OBJECT FOR M=2 

CALL SETUP 
CALL BUILD 

SET TIME LIMITS AND EVERY IP DATA POINT SAVED 

TSTART=0.0 

IP=IP*EXPFAC 

T=DT/2.+TSTART 

N=0 

CALL EADV 
DO 230 N=l, 192 
ADVANCE TIME 

T-T+DT/2. 

ADVANCE HFIELD 

CALL HADV 

ADVANCE TIME 

T=T+DT/2. 
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on no 


ADVANCE EFIELD 


C 

C 

C 


CALL EADV 
STORE FIELDS 

IF(MOD(N,IP) .EQ.O) CALL DATASAV 


230 CONTINUE 
C 

PRINT 250, T,N 

250 F ORMAT ( * 1 E XI T TIME(M=2)=*E12.3* .AFTER CYCLE*I4) 
CALL PRINOUT 
IF(HM.NE. 12) GO TO 260 
GO TO 500 
C 

260 CONTINUE 
C 

300 CONTINUE 
500 CONTINUE 
END 
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3.2 


GLOBAL VARIABLES 


The global variables for G3DXL are located in labeled common. The 
nine commons are labeled EFIELD, HFIELD, EXTRAS, GRID, UGRID, RAD, TSITEM, EBS, 
and OUT. The variables are described in Table 2. 

3.3 DATA FILES 

Input to Program G3DXL is made through data statements and by 
manually defining program variables, predominantly in subroutine SETUP. The 
variables input are identified in the discussion of the subroutines in which 
they are defined (Section 3.5). 

There is one output file for Program G3DXL. This printer output is 
referred to in the program as OUTPUT and is accessed via formated print 
statements. Descriptions of the output data are furnished in the discussions 
of the subroutines in which they are output (Section 3.5). 

3.H ERROR MESSAGES 

The only error message in G3DXL is in the main program DRIVER. The 
statement "MM INPUT ERROR" is printed if MM is not set to an allowed value. 
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Table 2. Global Variables 
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Table 2. (continued) 



31 



Table 2. (continued) 
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Table 2. (continued) 
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Table 2. (continued) 
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Table 2. (continued) 




>- 

M 

X 

M 

X 

>- 

i 

3 



X 

X 

>- 

>- 

PM 

PM 

£ £ 



OS 

os 

os 

OS 

OS 

OS 

*r- O 



*4- ^ 

<4- rT 

*4— 

4- 



CO *r— 
+-> 



O ^ 

O — - 

o ^ 

Ow 

o 

o — 

QJ u 









JC c i 



to 

CO 

to 

CO 

CO 

10 

3 



E c 

£ c 

E c 

B c 

£ c 

£ c 

M— 

CO 


C o 

C o 

C o 

C o 

C o 

C o 

S- 

c 

to 

(U *r— 

0) •*— 

QJ T- 

CU -1- 

CU •!“ 

CU *f- 

O CD 

o 

*a 

4-> +J 

4-> +J 

+-> -M 

4-> 

4-> -M 

+j +-> 


•f— 

c 

<0 

03 

03 

03 

03 

03 

•f— 

4-> 

o 

C 3 

C 3 

C 3 

C 3 

C 3 

C 3 

CO > 

Q. 

o 

•r- O" 

•*- O" 

*r- CD 

•r- cr 

*r- c r 

*r- cr 

0) *f- 

•f— 

CU 

LU 

LU 

LU 

LU 

LU 

LU 

+-> L 


CO 

"O 

X3 

TD 

*o 

XJ 

X7 

03 X3 

o 


(U o 

(U O 

<U O 

CU o 

<u o 

CU o 

C 

CO 


c +-> 

C +■* 

C -M 

C 4-> 

C -M 

c +-> 

•r- X3 

CU 

a> 

•r- 

•r* 

*r- 

•1— 

•r 

•f— 

U n* 

O 

3 

*4- CD 

»4- CD 

0- CD 

4- CD 

M- CD 

M- CD 

S- CU ^ 



a> c 

a> c 

<U C 

CU c 

<u c 

<u c 

O CO 


03 

-a *r- 

X3 -r- 

T3 

-o *- 

X3 -t- 

X3 *r- 

04- X 


> 

X3 

*o 

-a 

XJ 

X> 

-o 

(J C 



<u s- 

<U t- 

<U L- 

CU s- 

CU S- 

a> s- 

L- O 


(U 

r— O 

r- O 

r- O 

r- O 

*— o 

r— ■ O 

CU O U 


£ 

CD O 

CD U 

CD O 

CD U 

CD <J 

CD U 

E +-> cu 


•r— 

C O 

C O 

C CJ 

c a 

C tJ 

C CJ 

03 C0 


1— 

<C 03 

C 03 

C 03 

C 03 

<t 03 

C 03 

1- r- — 











CL CL 
PM PM 









*0 O 
NZZ 
O *— i 

ZLUX 








a> 

►— « •» •* 








c 

LU CL CL 








•r- CO 









CO <U 

>- o o 








Z C 

o z z 








•r— 

2S ►— « K— 4 








L -i-> 

• lu z 








CD 3 

1 1 1 ft ft > 








-C o 

** Q. O- «X 








4-> OS 

X X X GO 








o 

O O O C 

Q 

a 

Q 

o 

o 

o 

t— 


Z Z Z h- 

C 

c 


< 

<x. 

< 



►— « ►— « ►— < e£ 

OS 

os 

os 

OS 

OS 

OS 

Ll 


UJUJIQ 

LU 

LU 

LU 

LU 

LU 

Ui 

CD 

CD to 

£ 








c <u 

03 








c 

L- 








c •*— 

CD 

a_ 

CL 

a. 

Q- 

CL 

CL 

CL 

•«- 4-> 

C O 

z 

Z 

ID 

I D 

ID 

ID 

Z 

3 

■r- i- 

l— 

1 — 

H— 

h- 

h- 

h— 

h- 

aj o 

a3 a. 

LU 

LU 

LU 

LU 

LU 

LU 

LU 

q a: 

2= 

GO 

00 

CO 

CO 

GO 

GO 

GO 










CU 

r — 


* cm 

*» CNJ 

* CM 

" CM 

ft CM 

O 

tsi 


CO * 

CO • 

CO ** 

00 * 

CO •» 

CO ^ 

O 

•r— 


C\J CD 

CM CD 

CM CD 

CM CD 

CM CD 

CM CD 

CO 

GO 


CM 

CNJ 

- — CM 

CM 

" — CM 

— '<NJ 


C 

GO 








o 

< 







CL 


OS 







Z 

c 

\— 

Q 

Q 

O 

O 

CD 

o 

1— 

o 

X 

<C 


<c 

< 

< 

< 

LU 

o 

LU 

OS 

os 

OS 

OS 

OS 

OS 

GO 

cu 









r— 









-O 

03 CU 


>- 

hsl 

X 

fM 

X 

>- 

\— 

*r— £ 


X 

X 

>- 

>- 

PM 

r%i 

i 

&- ro 


z 

z 

Z 

z 

z 

z 

Ll 

<0 Z 
> 

1— 

h- 

1- 

I— 

1— 

h- 

i— 

1— 


35 




Table 2 . (continued) 
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3.5 


SUBROUTINES 


In this section the subroutines that comprise Program G3DXL are 
presented. A brief description of the function performed by the subroutine is 
provided along with a subroutine flow chart and a description of those 
variables local to the subroutine. 

SETUP - Problem space definition, initialization and parameter 
definitions 

BUILD - Scattering geometry definition 

EADV - Advances free space E-fields 

EBC - Sets electric field boundary conditions, lightning channels, 
and defines internal wires. 

ERAD - Outer radiation boundary condition 

HADV - Advances free space H-fields 
FIELDS - Defines indirect radiated fields 

SAVESB - Saves subboundary tangential E-fields for expansion technique 
OUTBND - Interpolates SAVESB fields for expanded calculation 
DATASAV - Saves data at selected test point locations 
PRINOUT - Outputs data from DATASAV, extends wire current data in time 
3.5.1 Subroutine SETUP 

The test configuration and conditions are specified in subroutine 
SETUP. The gridding for the cell space, the time step and the radiation 
factors are determined. The variables listed in Table 3 are set in the 
subroutine. A flow chart for Subroutine SETUP is furnished in Figure 3. Table 
H presents the variables local to the subroutine. 
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Table 3. Variables Input in Subroutine SETUP 


Variable 

Name 

Description 

DELX 

Length of cell where there is no cell expansion 
(x dimension) (m) 

DELY 

Width of cell where there is no cell expansion 
(z dimension) (m) 

DELZ 

Height of cell where there is no cell expansion 
(y dimension) (m) 

I DOWN 

I index of the first cell that is not expanded in 
the x direction 

IUP 

I index of the last cell that is not expanded in 
the x direction 

JDOWN 

J index of the first cell that is not expanded in 
the y direction 

JUP 

J index of the last cell that is not expanded in 
the y direction 

KDOWN 

K index of the first cell that is not expanded in 
the z direction 

KUP 

K index of the last cell that is not expanded in 
the z direction 

1ST 1 
JST ) 
KST J 

Starting point in the unexpanded grid (I,J,K) - 
( 1ST, JST, KST) , for the expanded grid, (I,J,K) = 
(0,0,0) 
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Figure 3* Flow Chart for Subroutine SETUP 
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Table 4. Local Variables - Subroutine SETUP 


Variable 

Name 

DTXI 

DTYI 

DTZI 

XPANX, XPANY 
XPANZ 

EXPFAC 

I DOWN 

IUP 

JDOWN 

JUP 

KDOWN 

KUP 

Q 

QXY 

QXZ 

QYX 

QYZ 

QZX 

QZY 

R1 

R2 



Minimum cell dimension in y 
Minimum cell dimension in z 


Cell expansion coefficient in x, y, z 
direction 

Expansion factor for M*2,4 loops 
See Table 3 


See Table 3 
See Table 3 


See Table 3 
See Table 3 


See Table 3 

Minimum value of QYZ, QZX, QXY, QZY, QXZ, and QYZ 

Minimum value of 6 's 

xy 

Minimum value of 9 ' s 

X z 

Minimum value of 9 1 s 

yx 

Minimum value of 9 ,'s 

yz 

Minimum value of 9 's 

Z X 

Minimum value of 9„, ,'s 

zy 

Temporary storage for distances in computation of 
radiation factors 
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The gridding of the cell space for the unexpanded loops is determined 
from DELX, DELY, and DELZ (see Table 3) and the expansion coefficients (XPANX, 
XPANY, and XPANZ) . The grid size within the nonexpanding area is simply DELX, 
DELY, and DELZ in the x, y, and z directions respectively. The majority, if 
not all, of the test item is usually contained within this area. Outside this 
area the cell size of each cell as one moves away from the center of the cell 
space is the previous cell size multiplied by the expansion factor for that 
coordinate, for example, 


DXOU+1) = DXO(I) * EXPANX 


The time step for the computations is defined as 


At 



(dy m ln ) ' 



-1 


where d& . is the minimum cell dimension for the i coordinate direction, 
mxn 

The radiation factors 0 ,0 ,0 ,0 ,0 ,0 ,R ,R ,R , 

yx zx xy zy xz yz yx zx xy 

R , R , and R (see Subroutine ERAD) are computed and stored for later used 
zy xz yz 

in the computations. 


Subroutine SETUP is called from DRIVER, the main program of G3DXL. 
Variables input and computed in Subroutine SETUP are printed before it returns 
to the program. 
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o n o o non non non non non 


SUBROUTINE SETUP 
C 

DIMENSION DDX0(210) ,DDY0(210) ,DDZ0(210) 

COMMON/EFIELD/EX(29 ,29 ,29) ,EY(29 ,29 ,29) ,EZ(29,29,29) 

COMMON/ HF IELD/ UX( 2 9 , 29 , 29 ) ,I1Y(29,29,29) ,HZ( 29, 29 , 29) 
COMMON/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ UGRID/ UX( 2 8 ) ,UY(28) ,UZ(28) ,UXO(29) ,UY0(29) ,UZ0(29) 

COMMON/ EXTRAS /NX, NY ,NZ , NX1 ,NY1 ,NZ1 , N,M,MQ,DT, XMU , EPSO.EPS .NPTLIM , 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, 

2 BETA, IDLS 

COMMON/RAD/RYX(28,29,2) ,RZX(29,28,2) ,RXY(28,29,2) ,RZY(29,28,2) , 

1 RXZ(28,29,2) ,RYZ(29,28,2) ,THYX(28,29,2) ,THZX(29,28,2) , 

2 THXY(28,29,2) ,THZY(29,28,2) ,THXZ(28,29,2) ,THYZ(29,28,2) 

COMMON/ EBS/EYXD( 28 ,29 , 3) ,EYXU(28 ,29 ,3) ,EZXD(29,28,3) ,EZXU(29,28, 3) 

1 ,EXYD(28,29,3) ,EXYU(28,29,3) ,EZYD(29,28,3) ,EZYU(29,28,3) , 

2 EXZD(28,29,3) ,EXZU(28,29,3) ,EYZD(29 ,28 , 3) ,EYZU(29,28,3) , 

3 N1,N2,N3 
C 

COMMON/CUMMUL/ CUMM( 24) ,L0C(24) 

C SET CYCLES BETWEEN PRINTS AND RECORD LENGTH 
C 

DATA LMAX/ 230/ 

SET EXPANSION FACTOR 
DATA EXPFAC/4.0/ 

SET PARAMETERS 

DATA PI, EPSO, XMU, C/3. 1415926536, 8. 854E-I2 , 1 . 256637E-6 , 3. E8/ 

SET PROBLEM SPACE DIMENSIONS 
DATA NX,NY,NZ ,NX1 ,NY1 ,NZ 1/29, 29, 29, 28, 28, 28/ 

SET AIRCRAFT SIGMA AND EPSILON 
DATA EPS, SIGMA/3. 1E-I1.1.0E2/ 

IDLS=0 FOR RADIATED FIELDS, IDLS=1 FOR DIRECT STRIKES 
DATA IDLS/ 1/ 

IP=2 

FIRST INITIALIZE FIELDS TO ZERO 

DO 10 1*1,29 
DO 10 J=l,29 
DO 10 K*1 ,29 
C 

EX(I,J,K)*0.0 

EY(I,J,K)*0.0 
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EZ(I,J,K)=0.0 
11X( I, J,K)=0.0 
11Y( I , J , K)=0 . 0 
HZ(I,J,K)=0.0 
10 CONTINUE 
C 

DO 20 1=1,28 
DO 20 J=1 , 29 
DO 20 K=1 , 3 
C 

EXYD( I , J ,K)=0 . 0 
EXYU ( I , J , K) =0 . 0 
EXZD(I,J,K)=0.0 
EXZU ( I , J , K.) =0 . 0 
EYXD( I, J ,K)=0. 0 
EYXU( I , J ,K)=0 . 0 
20 CONTINUE 
C 

DO 25 1=1,29 
DO 25 J=1 , 28 
DO 25 K=1 , 3 
EYZD( I , J ,K)=0. 0 
EYZU( I , J ,K.)=0 . 0 
EZXD( I, J,K)=0. 0 
EZXU(I,J,K)=0.0 
EZYD( I , J ,K)=0 . 0 
EZYU( I , J ,K)=0 . 0 
25 CONTINUE 
C 

N1=0 

N2=0 

N3=0 

C 

DO 30 LL=1 , 24 
CUMM(LL)=0.0 
30 CONTINUE 
C 

AMP=5 . 92E4 
ALPHA=4.08E6 
BETA=1 . 0E8 
C 

TX=0.002 

TY=0.005 

TZ=0.005 

C 

C 

C INPUT GRID FOR M-l 
C 

C SET THE UNEXPANDED CELL DIMEHSIOMS(DELX,DELY,DELZ) IN METERS, 

C EXPANSION FACTORS (XPANX.XPANY.XPANZ) AND THE LIMITS OVER 
C WHICH THE CELLS REMAIN CONSTANT(IUP,IDOWN, JUP, JDOWN,KUP,KDOWN) 
DELX=I . 2 
DELY=0.48 
DELZ=0.48 
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noon 


c 

XPANX=1 . 3 
XPANY=1 . 1 5 
XPANZ-1.15 
C 

IUP=25 

IDOWN-5 

JUP=25 

JD0WN=5 

K.UP-25 

KD0WN=5 

C 

DO 100 1=1 , IDOWN 
DX0( 1 ) =DELX* ( XPANX** ( IDOWN-I ) ) 
100 CONTINUE 
C 

DO 105 1= I DOWN, I UP 
DX0(I)=DELX 
105 CONTINUE 
C 

DO 110 I=IUP,NX1 
DX0(I)=DELX*(XPANX**(I+1-IUP) ) 
110 CONTINUE 
C 

DO 115 J=1 , JDOWN 
DYO ( J ) =DELY* ( XPANY** ( JDOWN- J ) ) 
115 CONTINUE 
C 

DO 120 J-JDOWN.JUP 
DY0( J)=DELY 
120 CONTINUE 
C 

DO 125 J=JUP,NY1 
DY0(J)=DELY*( XPANY** ( J+l-JUP) ) 
125 CONTINUE 
C 

DO 130 K= 1 , KDOWN 
DZO(K) =DELZ* (XPANZ** (KDOWN-K) ) 
130 CONTINUE 
C 

DO 135 K=KDOWN,KUP 
DZO(K.)=DELZ 
135 CONTINUE 
C 

DO 140 K=KUP,NZ1 
DZO(K) =DELZ* ( XPANZ** ( K+l -KUP ) ) 
140 CONTINUE 
C 

IF(M.EQ.l) GO TO 205 
INPUT GRID FOR M=2 


NX1=EXPF AC*NX1 
NY1=EXPFAC*NY1 
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NZ1=EXPFAC*NZ1 

NX=NX1+1 

NY=NY1+1 

Nz=NZl+l 

IST=7 

JST=9 

KST=12 

1 DOWN= ( EXPF AC* ( IDOWN- 1 ) )+l 
IUP=( EXPFAC*( IUP- 1 ) )+l 
XPANX-XPANX/ EXPFAC 

JDOWN= ( EXPFAC* ( JDOWN- 1 ) )+l 
JUP=( EXPFAC* ( JUP-1 ) )+l 
XPANY=XPANY/ EXPFAC 

KDOWN= ( EXPFAC* ( KDOWN-l ) )+l 
KUP=( EXPFAC* (KUP-1))+1 
XPANZ=XPANZ/ EXPFAC 

DELX=DELX/ EXPF AC 
DELY=DELY/EXPFAC 
DELZ=DELZ/ EXPFAC 

DO 145 1=1, IDOWN 
DDXO( I ) =DELX* ( XPANX** ( IDOWN-I ) ) 
145 CONTINUE 

DO 150 I=IDOWN, IUP 
DDXO( I)=DELX 
150 CONTINUE 

DO 155 I=IUP, NX1 
DDXO ( I ) =DELX* ( XPANX** ( I+l-IUP ) ) 
155 CONTINUE 

DO 160 J-l, JDOWN 

DDYO ( J ) =DELY* ( XPANY** ( JDOWN- J ) ) 

160 CONTINUE 

DO 165 J “JDOWN, JUP 
DDYO( J)=DELY 
165 CONTINUE 

DO 170 J=JUP,NY1 
DDYO ( J ) =DELY* ( XPANY** ( J+l-JUP ) ) 
170 CONTINUE 

DO 175 K=1 , KDOWN 
DDZO(K) =DELZ* ( XPANZ** ( KDOWN-K) ) 
175 CONTINUE 

DO 180 K=KDOWN,KUP 
DDZO(K)=DELZ 
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180 CONTINUE 


C 

DO 185 K=KUP,NZ1 
DDZO(K) =DELZ* ( XPANZ** ( K+l-KUP ) ) 
185 CONTINUE 
C 

IEXPST=( EXPFAC* ( 1ST- 1 ) )+l 
IEXPEND=IEXPST+27 
C 

DO 190 I=IEXPST , IEXPEND 
DXO( I-IEXPST+ 1 ) =DDXO ( I ) 

190 CONTINUE 
C 

JEXPST=(EXPFAC*( JST-1 ) )+l 
JEXPEND= JEXPST+2 7 
C 

DO 195 J= JEXPST , JEXPEND 
DY0(J-JEXPST+1)=DDY0(J) 

195 CONTINUE 
C 

KEXPST=(EXPFAC*(KST-1))+1 

KEXPEND=KEXPST+27 

C 

DO 200 K=KEXPST , KEXPEND 
DZ0(K-K£XPST+1)=DDZ0(K) 

200 CONTINUE 
C 

NX1=NX1/EXPFAC 
NY 1=NY1/ EXPFAC 
NZ1=NZ1/ EXPFAC 
NX=NX1+1 
NY=NY1+1 
NZ=NZ1+1 
C 

ID0WN=((ID0WN-1)/EXPFAC)+1 
IUP=((IUP-1)/EXPFAC)+1 
JDOWN=( (JDOWN— 1) /EXPFAC)+1 
JUP=( ( JUP— 1)/EXPFAC)+1 
KDOWN=( (KDOWN-1 ) /EXPFAC)+1 
KUP=( (KUP-1)/ EXPFAC )+l 
C 

205 CONTINUE 
C 

D=0. 0 

DO 210 1=1 ,NX1 
D=D+DXO( I ) 

210 CONTINUE 
C 

DO 215 1=2, NX 
X0( l)=-(D/2 . ) 

IF(M.EQ.2) X0( 1)=UX0( 1ST) 
X0(I)=X0(I-1)+DX0(I-1) 
IF(M.EQ.2) GO TO 215 
UX0(1)=X0(1) 

UXO(I)=XO(I) 


46 



215 CONTINUE 


C 

W=0.0 

DO 220 J=1,NY1 
IJ=W+DY0(J) 

220 CONTINUE 
C 

DO 225 J=2 ,NY 
Y0(l)=-(W/2.) 

IF(M.EQ. 2) Y0( 1)=UY0(JST) 
Y0( J)=YO( J— 1)+DY0( J— 1 ) 
IF(M.EQ. 2) GO TO 225 
UYO( 1)=Y0( 1) 

UYO(J)=YO(J) 

225 CONTINUE 
C 

H=0 . 0 

DO 230 K=1,NZ1 
H=H+DZO(K) 

230 CONTINUE 
C 

DO 235 K=2 , NZ 
Z0(l)=-(H/2.) 

IF(M.EQ. 2) Z0( 1)=UZ0(KST) 
Z0(K)=Z0(K-1)+DZ0(K-1) 
IF(M.EQ.2) GO TO 235 
UZ0(1)=Z0(1) 

UZO(K)=ZO(K) 

235 CONTINUE 
C 

DO 240 1=1, NX1 
X(I)=(X0(I)+X0(I+l))/2. 
IF(M.EQ. 2) GO TO 240 
UX(I)=X(I) 

240 CONTINUE 
C 

DO 245 J=1 ,NY1 
Y(J)=(Y0(J)+Y0(J+l))/2. 
IF(M.EQ. 2) GO TO 245 
UY( J)=Y(J) 

245 CONTINUE 
C 

DO 250 K=1 ,NZ1 
Z(K)=(Z0(K)+Z0(K+l))/2. 
IF(M.EQ. 2) GO TO 250 
UZ(K)=Z(K) 

250 CONTINUE 
C 

DO 260 1=2, NX1 
DX(I)=X(I)-X(I-1) 

260 CONTINUE 

DX( 1)=DX(2)*XPANX 
DX(29)=DX(28)*XPANX 
C 

DO 265 J=2 ,NY1 
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DY(J)=Y(J)-Y(J-1) 

265 CONTINUE 

DY ( 1 ) =DY ( 2 ) *XPANY 
DY(29)=DY(28)*XPANY 
C 

DO 270 K=2 ,NZ1 
DZ(K)=Z(K)-Z(K-1) 

270 CONTINUE 

DZ ( 1 ) =DZ ( 2 ) *XPANZ 
DZ(29)=DZ(28)*XPANZ 
C 

DO 275 1=1, NX1 
DXOI(1)=1./DXO(I) 

275 CONTINUE 
C 

DO 280 J=1,NYI 
DY0I(J)=1./DY0(J) 

280 CONTINUE 
C 

DO 285 K=1,NZ1 
DZ0I(K)=1 ./DZO(K) 

285 CONTINUE 
C 

DO 290 1=1, NX 
DXI ( I ) = I . / DX ( I ) 

290 CONTINUE 
C 

DO 295 J=1 ,NY 
DYI(J)=1./DY(J) 

295 CONTINUE 
C 

DO 300 K=1 ,NZ 
DZI(K)=1./DZ(K) 

300 CONTINUE 
C 

DTXI=C/DX0(NXl/2) 

DTYI=C/DY0(NYl/2) 

DTZI»C/DZ0(NZl/2) 

DT=1 . /SQRT(DTXI**2+DTYI**2+DTZI**2) 

C 

IF(M.EQ. 2) GO TO 199 

EYX RADIATION FACTORS 
OYX = 1. 

DO 61 K=1,NZ 
DO 61 J-l.NYl 
XI = X0( 1)**2 

Y2 = Y( J)**2 

Z2 = Z0(K)**2 

R1 = SQRT(X1 + Y2 + Z2) 

X2 = X0(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 

T11YX( J ,K, 1) = 1. - (R1-R2)/ (C*DT) 

IF(THYX(J,K,1).LT.QYX) QYX = THYX(J,K,1) 
RYX( J ,K, 1) = R2/R1 
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XI = X0(NX)**2 
R1 = SQRT(X1 + Y2 + Z2) 

X2 = X0(NX1)**2 
R2 = SQRT(X2 + Y2 + Z2) 

THYX(J,K,2) = 1. - (R1-R2)/ (C*DT) 

IF(THYX(J,K,2) .LT.QYX) QYX = THYX(J,K,2) 
RYX(J,K,2) = R2/R1 
bl CONTINUE 

EZX RADIATION FACTORS 

QZX = 1. 

DO 62 K=1 , NZ1 

DO o2 J=L ,NY 
XI = XO( 1)**2 

Y2 = YO( J)**2 

Z2 = Z(K)**2 

R1 = SQRT(X1 + Y2 + Z2) 

X2 = X0(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THZX(J,K,1) = 1. - (R1-R2)/ (C*DT) 

IF(THZX( J,K, 1) .LT.QZX) QZX = THZX(J,K,1) 
RZX(J,K,1) = R2/R1 
XI = XO(NX)**2 
R1 = SQRT(X1 + Y2 + Z2) 

X2 = XO(NXl)**2 
R2 = SQRT(X2 + Y2 + Z2) 

THZX(J,K,2) = 1. - (R1-R2)/ (C*DT) 

IF(TIiZX(J,K, 2). LT.QZX) QZX = THZX(J,K,2) 
RZX(J,K,2) - R2/R1 

62 CONTINUE 

EXY RADIATION FACTORS 

QXY = 1. 

DO 63 K=1 ,NZ 

DO 63 1=1, NX1 

X2 = X(I)**2 

Y1 = YO(l)**2 

Z2 = ZO(K)**2 

R1 = SQRT(X2 + Y1 + Z2) 

Y2 = YO(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THXY ( I , K , 1 ) = 1. - (R1-R2)/ (C*DT) 

IF(TilXY( I , K, 1) .LT.QXY) QXY = THXY(I,K,1) 
RXY ( I , K , 1 ) = R2/R1 
Y1 = YO(NY)**2 
R1 = SQRT(X2 + Y1 + Z2) 

Y2 = YO(NYl)**2 
R2 = SQRT(X2 + Y2 + Z2) 

THXY ( I , K , 2 ) = 1. - (R1-R2)/ (C*DT) 

IF(THXY(I,K, 2). LT.QXY) QXY = THXY(I,K,2) 
RXY( I,K, 2) = R2/R1 

63 CONTINUE 
C 
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EZY RADIATION FACTORS 


QZY = I. 

DO 64 K=l,NZl 
DO 64 1=1, NX 

X2 = X0(I)**2 

Y1 = Y0(l)**2 

Z2 = Z(K)**2 

R1 = SQRT(X2 + Y1 + Z2) 

Y2 = Y0(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THZY( I ,K, 1) = I. - (RI-R2)/ (C*DT) 
IF(THZY(I,K, 1) .LT.QZY) QZY = THZY(I.K.l) 
RZY(I,K, 1) = R2/R1 
Y1 = Y0(NY)**2 
R1 = SQRT(X2 + Y1 + Z2) 

Y2 = Y0(NY1)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THZY( I ,K, 2) = 1. - (R1-R2)/ (C*DT) 

IF(THZY(I,K, 2) .LT.QZY) QYZ=THYZ(I,K,2) 
RZY( I,K,2)=R2/R1 
64 CONTINUE 

EXZ RADIATION FACTORS 


QXZ = I. 

DO 65 1=1, NX1 

DO 65 J=1 ,NY 

X2 = X(I)**2 

Y2 = Y0( J)**2 

Z1 = Z0(l)**2 

R1 = SQRT(X2 + Y2 + Zl) 

Z2 = Z0(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THXZ ( I , J , 1 ) = 1. - (R1-R2)/ (C*DT) 

IF(THXZ(I,J,1) .LT.QXZ) QXZ = THXZ(I,J,1) 
RXZ ( I , J , 1 ) = R2/R1 
Zl = Z0(NZ)**2 
R1 = SQRT(X2 + Y2 + Zl) 

Z2 = Z0(NZ1)**2 

R2 = SQRT(X2 + Y2 + Z2) 

THXZ( I , J , 2) = 1. - ( R1-R2 ) / ( C*DT ) 

IF(TIiXZ(I,J, 2). LT.QXZ) QXZ = THXZ(I,J,2) 
RXZ ( I , J , 2 ) = R2/R1 
65 CONTINUE 

EYZ RADIATION FACTORS 


QYZ = 1. 

DO 66 J=1 ,NY1 

DO 66 1=1, NX 

X2 = X0(I)**2 

Y2 = Y(J)**2 

Zl = Z0( 1)**2 

R1 = SQRT(X2 + Y2 + Zl) 

Z2 = Z0(2)**2 

R2 = SQRT(X2 + Y2 + Z2) 
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THYZ(I,J,l)=l.-(Rl-R2)/(C*DT) 

IF(THYZ(I,J,1).LT.QYZ) QYZ“THYZ( I , J, 1 ) 

RYZ(I, J,1)=‘R2/R1 

Z1»Z0(NZ)**2 

R1=SQRT(X2+Y2+Z1) 

Z2=Z0(NZ1)**2 

R2=*SQRT(X2+Y2+Z2) 

THYZ(I,J,2) - 1. - (R1-R2)/ (C*DT) 

IF(THYZ(I,J,2) .LT.QYZ) QYZ - THYZ(I,J,2) 
RYZ(I,J,2) = R2/R1 
66 CONTINUE 

Q = AMIN1(QYX,QZX ,QXY,QZY , QXZ , QYZ) 

MQ = -Q/2. 

PRINT VARIABLES 


199 PRINT 201 

PRINT 202, DT, EPS, SIGMA 
PRINT 99 

PRINT 104, NX, NX1 , NY, NY1, NZ, NZ1 
PRINT 101 

PRINT 98, (XO(I), I = 1 ,NX) 

PRINT 102 

PRINT 98, (Y0( J) , J - 1 ,NY) 

PRINT 107 

PRINT 98, (DXO(I), I - 1,NX1) 

PRINT 103 

PRINT 98, ( ZO(K) , K = 1,NZ) 

PRINT 106 

PRINT 98, QYX,QZX,QXY,QZY,QXZ,QYZ,Q 


FORMAT STATEMENTS 


98 FORMAT (10F10.4) 

101 FORMAT (*0X0 GRID*) 

102 FORMAT (*0Y0 GRID*) 

103 FORMAT (*0Z0 GRID*) 

104 FORMAT (2015) 

99 FORMAT ( 5110 NX,2X,*NXl NY NY I NZ NZI MQ*) 

106 FORMAT ( 1H0 , 5X , *QYX* , 6X , *QZX* , 6X , *QXY* , 6X , *QZY* , 6X , *QXZ* , 6X , *QYZ* , 

1 7X,*Q*) 

107 FORMAT (*ODXO GRID*) 

202 FORMAT (3EI2.4) 

201 FORMAT (*1DT, EPS, SIGMA*) 

390 F0RMAT(2F10.3,2EI2.4) 

400 FORMAT (3F 10.0) 

440 FORMAT(6I5) 

RETURN 

END 
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3.5.2 Subroutine BUILD 

The NOPE array, which defines the test time, is filled in Subroutine 
BUILD. Each cell in the cell space is identified by three indices (I,J,K) (see 
Figure 1). The entire array is initialized to zero. For perfect conductors, 
if the I,J,K^ cell is a three-dimensional cell in the test time, NOPE(I,J,K) 
is set equal to four (4). If the YZ plane adjacent to vertex I,J,K of the cell 
is a two-dimensional piece of the test time, NOPE(I,J,K) is set equal to one 
(1). For the XZ plane, NOPE(I,J,K) is set to two (2) and to three (3) if the 
XY plane is a two-dimensional part of the test item. Where planes meet NOPE is 
set to 12 (YZ and XZ planes), 13 (YZ and XY planes, 23 (XZ and XY planes) and 
123 (YZ, XZ and XY planes). When an imperfect conductor is modeled 4+8 and 

1.2.3 5,6,7. The values of the NOPE array for each XZ plane cut of the test 
item are printed and can be used to check the test item structure. This 
routine is replaced in its entirety whenever a computation on a new structure 
is required. 

The version of Subroutine BUILD included here builds the F-106B 
aircraft. A flow chart for the routine is given in Figure 4. Subroutine BUILD 
is called by the main program DRIVER. 
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SUBROUTINE BUILD 


C 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 

1 NN,NPTS,LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
COMMON/TSITEM/NOPE(29 ,29,29) 

C 

DO 100 1=1, NX 
DO 100 J=1 ,NY 
DO 100 K=1,NZ 
NOPE( I, J,K)=0 
100 CONTINUE 

CYLINDER TEST CASE FOR CODE OPERATION CHECK 
TEST=0 

IF(TEST.EQ.O) GO TO 8888 

BUILD M=1 GEOMETRY 
IF(M.NE.l) GO TO 200 
DO 110 1=5,24 
DO 110 J=13,16 
DO 110 K=13, 16 
N0PE(I,J,K)-8 
110 CONTINUE 
GO TO 700 

1100 CONTINUE 

DO 1101 1=5,25,20 
DO 1101 J=13 , 16 
DO 1101 K=13,16 
N0PE( I, J ,K)=5 

1101 CONTINUE 

DO 1102 1=5,24 
DO 1102 J=13,17,4 
DO 1102 K=13 , 16 

IF(NOPE(I,J,K).EQ.l) GO TO 1111 
NOPE( I, J,K)=6 
GO TO 1102 

1111 N0PE(I, J,K)=56 

1102 CONTINUE 

DO 1103 1=5,24 
DO 1103 J=13,16 
DO 1103 K=13 , 17 , 4 
IF(NOPE(I,J,K) .EQ.l) GO TO 1112 
IF(NOPE( I , J ,K) .EQ.2) GO TO 1113 
IF(NOPE(I,J,K).EQ.l 2) GO TO 1114 
N0PE(I, J,K)=7 
GO TO 1103 

1112 N0PE( I, J,K)=57 
GO TO 1103 

1113 NOPE(I,J,K)=67 
GO TO 1103 

1114 NOPE(I,J,K)=567 

1103 CONTINUE 
GO TO 700 
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C BUILD M=2 GEOMETRY 
200 CONTINUE 

IF(M.NE.2) GO TO 300 
DO 120 1=1,28 
DO 120 J=5,20 
DO 120 K=5,20 
N0PE( I , J ,K)=8 
120 CONTINUE 
C 

GO TO 700 
C 

8888 CONTINUE 
C 

BUILD M=1 F 106 B GEOMETRY 

IF(M.NE.l) GO TO 1200 

FUSELAGE 

DO 10 1=8,20 
DO 10 J=ll , 12 
DO 10 K=13, 16 
N0PE( I , J ,K)=4 

10 CONTINUE 

DO 11 1=11,22 
DO 11 J=13 , 14 
DO 11 K=13, 16 
N0PE(I, J,K)=4 

11 CONTINUE 

N0PE(9 , 13 , 14)=4 
N0PE(9 , 13 , 15)=4 

NOPE (10, 13, 13) =4 
NOPE ( 10 , 13 , 14)=4 
NOPE(10,13,15)=4 
NOPE (10, 13, 16) =4 

N0PE(10,14,14)=4 
NOPE (10, 14, 15) =4 

DO 12 1=21,22 
J=12 

DO 12 K=13 , 16 
N0PE( I , J ,K)=4 

12 CONTINUE 

DO 13 1=13,15 
DO 13 J=12 , 13 
DO 13 K=12 , 17,5 
NOPE( I, J ,K)=4 

13 CONTINUE 

WINGS 
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c 

DO 14 1=14,20 
J=12 

DO 14 K=9 ,20 
N0PE( I , J ,K)=4 

14 CONTINUE 
C 

NOPE ( 14 , 12 , 10)=0 
HOPE (14,12,19)=0 
NOPE (14, 12, 9 )=0 
N0PE(14, 12 ,20)=0 
NOPE(15,12,9)=0 
NOPE ( 1 5, 12,20) =0 
C 

DO 15 1=18,20 
J=12 

DO 15 K=7 ,22 
N0PE( I , J ,K)=4 

15 CONTINUE 
C 

N0PE(18,12,7)=0 
NOPE (18, 12, 22) =0 

STABILIZER 

DO 16 1=19,20 
DO 16 J-15,17 
K=15 

N0PE( I, J ,K)=3 

16 CONTINUE 

N0PE(19,17,15)=0 

DO 17 1=21,22 
DO 17 J=15,21 
K=15 

N0PE(I,J,K)=3 

17 CONTINUE 

N0PE(21 ,20, 15)=0 
H0PE(21 ,21 , 15)=0 

GO TO 700 

1200 CONTINUE 

BUILD M=2 GEOMETRY 

IF(M.NE. 2) GO TO 1300 
C 

C BUILD SOLID PORTION FIRST 
C 

C TREAT SOLID PORTION ON SUB BOUNDARY 
C 

1=28 
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DO 20 J=9,24 
DO 20 K=5,20 
NOPE( I, J,K)=4 

20 CONTINUE 
C 

DO 21 J=I3,20 
DO 21 K=l,24 
N0PE(I,J,K)=4 

21 CONTINUE 

TREAT SOLID UPPER BODY 

DO 22 1=17,27 
DO 22 J=15,24 
DO 22 K=6 ,19 
N0PE(I,J,K)=4 

22 CONTINUE 
C 

DO 23 1=13,16 
DO 23 J=13,20 
DO 23 K=6,19 
NOPE( I, J,K)=4 

23 CONTINUE 
C 

DO 24 1=13,16 
DO 24 J-21,24 
DO 24 K=7 ,18 
N0PE(I,J,K)=4 

24 CONTINUE 
C 

DO 25 1=10,12 
DO 25 J=13, 18 
DO 25 K= 6 , 1 9 
NOPE( I, J ,K)=4 

25 CONTINUE 
C 

1=12 

DO 26 J=19 ,23 
DO 26 K=8 , 1 7 
N0PE(I,J,K)=4 

26 CONTINUE 


DO 27 J=19,22 
DO 27 K=9 ,16 
NOPE( I , J ,K)=4 

27 CONTINUE 
C 

1=10 

DO 28 J=19 ,20 
DO 28 K=10, 15 
NOPE(I, J ,IC)=4 

28 CONTINUE 
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DO 29 J=18,19 
DO 29 K=ll , 14 
N0PE(I,J,K)=4 

29 CONTINUE 
C 

DO 30 1=6,9 
J=17 

DO 30 K=7 , 18 
H0PE(I,J ,K)=4 

30 CONTINUE 
C 

DO 31 1=4,5 
J=16 

DO 31 K=8,17 
N0PE( I, J ,K)=4 

31 CONTINUE 
C 

C TREAT Y-Z THIN PLANES(NOPE(I,J,K)=l) 

C 

C BULKHEAD B1 

C 

1=4 

DO 32 J=9 ,15 
DO 32 K=8, 17 
NOPE(I,J,K)=l 

32 CONTINUE 

HOLES IN B1 

DO 33 J=9,14 
DO 33 K=12,13 
N0PE(I,J,K)=0 

33 CONTINUE 

BULKHEAD B2 


1=6 

DO 34 J=9 , 16 
DO 34 K=7 ,18 
N0PE(I, J,K)=1 

34 CONTINUE 

HOLES IN B2 

DO 35 J=9, 10 
DO 35 K=14 , 17 
NOPE(I, J ,K)=0 

35 CONTINUE 

BULKHEAD B4 


1=30 

DO 36 J=9 , 12 
DO 36 K=6 ,19 
N0PE( I, J,K)=1 
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36 CONTINUE 

HOLES IN B4 

DO 37 J=9 , 10 
DO 37 K=17 , 18 
N0PE( I , J ,K)=0 

37 CONTINUE 

BULKHEAD B6 
1=16 

DO 38 J=9 , 14 
DO 38 K=6 , 19 
N0PE( I, J ,K)=1 

38 CONTINUE 

HOLES IN B6 

DO 39 J=9 , 10 
DO 39 K=18 , 19 
N0PE(I,J,K)=0 

39 CONTINUE 

BULKHEAD B8 
1=20 

DO 40 J=9 , 14 
DO 40 K=6, 19 
N0PE(I,J,K)=1 

40 CONTINUE 

HOLES IN B8 

DO 41 J=13 , 14 
DO 41 K=18 , 19 
N0PE( I , J ,K)=0 

41 CONTINUE 

TREAT X-Z THIN PLANES (N0PE( I, J,K)=2) 

DO 42 1=4,27 
J=9 

DO 42 K=6 , 19 

IF ( NOPE ( I , J , K) . EQ . 1 ) GO TO 420 
N0PE( I , J ,K)=2 
GO TO 42 

420 NOPE( I , J ,K)=12 

42 CONTINUE 
C 

DO 43 1=4,9 
DO 43 K=6,19,13 
N0PE(I,J,K)=0 

43 CONTINUE 
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DO 44 1=4,5 
DO 44 K=7 ,18,11 
NOPE( I, J ,K)=0 

44 CONTINUE 

TREAT X-Y THIN PLANES (NOPE( I ,J ,K)=3) 

DO 45 1=10,15 
DO 45 J=9,12 
DO 45 K=6 ,20,14 
IF(NOPE(I,J,K) .EQ.l) GO TO 450 
IF(NOPE(I,J,K).EQ.2) GO TO 451 
IF(N0PE(I,J,K).EQ.12) GO TO 452 
NOPE( I, J ,K)=3 
GO TO 45 

450 NOPE( I, J ,K)=13 
GO TO 45 

451 NOPE(I,J,K)=23 
GO TO 45 

452 NOPE(I, J,K)=123 

45 CONTINUE 

DO 46 1=16,27 
DO 46 J=9 , 14 
DO 46 K=6,20,14 
IF(NOPE(I, J,K) .EQ.l) GO TO 460 
IF(NOPE(I,J,K) .EQ.2) GO TO 461 
IF(N0PE(I,J,K).EQ.12) GO TO 462 
NOPE( I, J ,K)=3 
GO TO 46 

460 N0PE(I,J,K)=13 
GO TO 46 

461 NOPE(I,J,K)=23 
GO TO 46 

462 NOPE( I, J,K)=123 

46 CONTINUE 

DO 47 1=6,9 
DO 47 J=9 , 16 
DO 47 K=7 ,19,12 
IF(NOPE(I,J,K) .EQ.l) GO TO 470 
IF(N0PE(I,J,K) .EQ.2) GO TO 471 
IF(NOPE(I,J,K) .EQ.12) GO TO 472 
N0PE(I,J,K)=3 
GO TO 47 

470 NOPE( I, J ,K)=13 
GO TO 47 

471 NOPE( I, J ,K)=23 
GO TO 47 

472 NOPE( I , J ,K)=123 

47 CONTINUE 

DO 48 1=4,5 
DO 48 J=9,15 
DO 48 K=8 ,18,10 
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IF(NOPE(I,J,K).EQ.l) GO TO 480 
IF(N0PE(I,J,K).EQ.2) GO TO 481 
IF(NOPE(I,J,K) .EQ.12) GO TO 482 
N0PE(I,J,K)=3 
GO TO 48 

480 NOPE( I, J ,K)=13 
GO TO 48 

481 NOPE(I,J,K)=23 
GO TO 48 

482 HOPE(I,J,K)=123 
48 CONTINUE 

C 

GO TO 700 
C 

401 F0RMAT(1H1 , / / ,5X,7HNOPE(I ,12 ,3U, K) ,// ) 

402 FORMAT (2 014) 

700 CONTINUE 

RETURN 

END 


61 


3 . 5.3 


Subroutine EADV 


In subroutine EADV, the equation 


-* -> 9E -► 

VxH=e~+aE 


( 1 ) 


-» -* ■+ 

is applied to the scattered fields and is used to find E from E and H at 
earlier times. This "time stepping" solution is found for each component 
of E, i.e., Equation (1) is written explicitly for each component as 


3E 

— - + 

oE 

3H 
2 

3H 

_I 

3t 

X 

3y 

3z 

3E 

_1 + 

cjE 

* x | 

li 
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y 

3z 

3x 

3E 

— — + 

oE 

3H 

. _J£ _ 
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X 

3t 

0 z 

3x 
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These equations are cast in their finite-difference form and solved for E at 
the latest time. Thus 


X X 


( 2 ) 


+ (At/e) 


H"(I,J f K) - H"(I,J-1 ,K) 
z z 

Y(J) - Y(J-1) 


Hy(I,J,K) - H"(I,J,K-1) 
Z(K) - Z(K-1) 


_n+1/2._ _ _n-1/2,_ _ ... . . 

Ey (I ,J ,K) ** Ey (I , J , K) + (At/e) 


( 3 ) 


/ h"(i,j, 

( — 


K) - H^(I , J ,K~1 ) 
Z (K) - Z(K-1) 


H"(I f J,K) - H"(I-1 ,J,K) \ 
X(I) - X(I-1) J 
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E n+1/2 (I,J,K) = E n “ 1/2 (I,J,K) + (At/e) 
z z 


( 4 ) 


Hy(I,J,K) -h"(I-1,J,K) h"(I,J,K) - Hjf(I,J-1,K) 

X(I)-X(I-1) Y(J) ~ Y(j-1) 

EADV is called every program cycle by the main program, aa is the 
Subroutine HADV for advancing the H-fielda. It uaes the prior valuea of the 
fields to calculate the new values. Initially, these prior values are zero. 

The first nonzero valuea are thoae valuea for E and H that imposed at the 
scatterer's boundary. These values are accessed by a call to EBC after the 
Equations (2), (3), (4) have been advanced in EADV. Thus, on the next program 
cycle these values will be available and nonzero values of the fields off the 
boundary will start appearing when Equations (2), (3), and (4) are advanced 
again. Since EBC is always called after the fields are advanced, the proper 
boundary conditions are maintained. 

Additionally, EADV stores the E-field componenta one cell in from 
each face of the outer boundary that are used in the calculation of the E-field 
component values on the outer boundary. This calculation, based on the 
radiation boundary condition [2], is performed in Subroutine ERAD. 
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Figure 5. Flow Chart for Subroutine EADV 





non non non 


SUBROUTINE EADV 


COMMON/EFIELD/EX(29 , 29 , 29) ,EY(29,29,29) ,EZ(29,29,29) 
COMMON/HFIELD/HX(29,29,29) ,HY (29, 29, 29), HZ (29, 29,29) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , N Y 1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 

1 NN, NPTS.LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,X0(29) , Y0(29) , Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,LY0I(28) ,UZ0I(28) 

COMMON/ EBS/EYXD(28 ,29,3) ,EYXU(28,29,3) ,EZXD(29,28,3) ,EZXU(29,28,3) 

1 ,EXYD(28,29,3) ,EXYU(28,29 ,3) ,EZYD(29,28,3) ,EZYU(29,28,3) , 

2 EXZD(28,29,3) ,EXZU(28,29,3) ,EYZD(29,28,3) ,EYZU(29,28,3) , 

3 N1,N2,N3 

COMMON/TSITEM/NOPE( 29 ,29,29) 

C 

DTE=DT/EPSO 

ADVANCE EX 

DO 1 I = 1 , NX1 
DO 1 J = 2 ,NYI 
DO 1 K = 2 ,NZ1 

IF(NOPE(I, J,K) .NE.O.AND.NOPE(I,J,K).NE.1.AND.NOPE(I,J,K).NE.5) 

1 GO TO 1 

EX( I , J ,K)=EX( I , J , K)+DTE* ( (HZ( I , J ,K)-HZ( I , J-l , K) )*DYI( J) 

1 -(HY(I,J,K)-HY(I,J,K-1))*DZI(K)) 

1 CONTINUE 

ADVANCE EY 

DO 2 I = 2 ,NX1 
DO 2 J = 1 ,NY1 
DO 2 K = 2 ,NZ I 

IF(NOPE(I, J,K) . NE . 0 . AND . NOPE (I,J,K).NE.2. AND. NOPE( I,J,K).NE.6) 

1 GO TO 2 

EY(I,J,K)=EY(I,J,K)+DTE*((HX(I,J,K)-liX(I,J,K-l))*DZI(K) 

1 -(HZ( I , J ,K)-IIZ( 1-1 , J ,K) )*DXI( I) ) 

2 CONTINUE 

ADVANCE EZ 

DO 3 I = 2 , NXI 
DO 3 J = 2 ,NY1 
DO 3 K = 1 ,NZ1 

IF(NOPE(I, J,K) .NE.O.AND.NOPE(I,J,K).NE.3.AND.NOPE(I,J,K).NE.7) 

1 GO TO 3 

EZ(I,J,K)=EZ(I,J,K)+DTE*((HY(I,J,K)-HY(I-L,J,K))*DXI(I) 

1 -(HX( I, J,K)-HX( I, J-l ,K) )*DYI( J) ) 

3 CONTINUE 

IF (M.EQ.2) GO TO 10 
IF(N2.NE.N3) CALL ERAD 
10 CONTINUE 
CALL EBC 

IF (M.EQ.2) CALL OUTBND 
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IF (M.EQ.2) RETURN 
IF (MOD(N,MQ) .NE.O) RETURN 
N3=N2 
N2=N1 
N1=N 
C 
C 

C ADVANCE EXY 
C 

DO 11 I = 1 ,NXI 
DO 11 K = 2 ,NZ1 
EXYD(I,K,3) = EXYD( I , K, 2 ) 
EXYD(I,K,2) = EXYD( I,K, 1) 
EXYD(I,K,1) = EX(I,2,K) 
EXYU(I,K, 3) = EXYU(I,K,2) 
EXYU(I,K,2) = EXYU(I,K,I) 
EXYU(I,K,I) = EX(I,NY1,K) 

II CONTINUE 
C 

C ADVANCE EXZ 
C 

DO 111 I = 1 , NX1 
DO 111 J = 2 ,NY1 
EXZL(I,J,3) = EXZD(I,J,2) 
EXZD( I , J , 2) = EXZD(I,J,1) 
EXZD(I,J,1) = EX(I,J,2) 

EXZU( I , J , 3) = EXZU(I,J,2) 
EXZU(I,J,2) = EXZU(I,J,1) 
EXZU(I, J, 1) = EX( I , J ,NZ1 ) 

111 CONTINUE 
C 

C ADVANCE EYX 
C 

DO 22 J = 1,NY1 

DO 22 K = 2 ,NZ1 

EYXD( J ,K, 3) = EYXD(J,K,2) 

EYXD(J,K,2) = EYXD(J,K, 1) 

EYXD(J,K, 1) = EY(2,J,K) 

EYXU( J,K, 3) = EYXU(J,K,2) 
EYXU(J,K,2) = EYXU(J,K, 1) 
EYXU( J , K, 1) = EY( NX1 , J ,K) 

22 CONTINUE 
C 

C ADVANCE EYZ 
C 

DO 222 I = 2,NX1 
DO 222 J = 1,NY1 
EYZD(I,J,3) = EYZD(I,J,2) 
EYZD( 1 , J , 2) = EYZD( I, J, 1) 
EYZD(I, J , 1) = EY(I, J, 2) 
EYZU(I,J,3) » EYZU( I, J,2) 
EYZU(I,J,2) = EYZU( I, J , 1) 
EYZU( I , J, 1) = EY(I,J,NZ1) 

222 CONTINUE 
C 
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ADVANCE EZX 

DO 33 J - 2 , NY 1 
DO 33 K - 1 ,NZ1 
EZXD(J,K,3) - EZXD(J,K,2) 
EZXD(J,K,2) - EZXD( J ,K, 1 ) 
EZXD(J,K,1) - EZ(2,J,K) 
EZXU( J ,K, 3) = EZXU(J,K,2) 
EZXU( J , K, 2) - EZXU(J,K,1) 
EZXU(J,K, 1) = EZ(NX1,J,K) 
33 CONTINUE 


ADVANCE EZY 


C 


DO 333 I = 
DO 333 K = 
EZYD(I ,K, 3) 
EZYD( I , K, 2) 
EZYD( I , K, 1 ) 
EZYU(I,K,3) 
EZYU(I,K,2) 
EZYU( I , K, 1 ) 
333 CONTINUE 


,NX1 

,NZ1 

= EZYD(I,K,2) 
=* EZYD( I ,K, 1) 
- EZ(I,2,K) 

= EZYU(I,K,2) 
= EZYU(I,K, 1) 
= EZ( I,NY1 ,K) 


RETURN 

END 
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Subroutine EBC 


Subroutine EBC sets the appropriate electric field boundary 
conditions for the test item. On the surface of a perfectly conducting 
aircraft, the total tangential E-field vanishes so 

gtan _ _gtan 

scattered incident 

The conditions are set in one part of the subroutine for any cell that defines 

a three-dimensional portion of the test item - E , E and then E . The setting 

x y z 

of boundary conditions for the two dimensional surfaces is done in a separate 
part of the subroutine. 


The values of the NOPE array for the cell and its surrounding cells 

are used to determine what the E-field value should be. If NOPE(I,J,K) is 

zero, no action is taken. For a cell with NOPE(I,J,K) = 4, the three 

surrounding cells are checked. If any one of the three is missing, the field 

tan 

is on the surface and E is set to -E. . . If all three cells are present E 

incident 

is not set, since it is then internal to the test item. Figure 6 indicates the 

situation for E x on cell I,J,K. It is clear that if either cell ©© °r ©is 

not a part of the test item, E x for cell© would be on an external surface of 

the item. Similar situations can be defined for E and E . 

y z 

In addition, if, for E x , cell I,J,K is the outermost cell in either 

the J or the K directions, the E values at the far edges of this last cell are 

t an ^ 

set equal to — E inQidenfc . The same is done for E y , if the cell is the outermost 

in either the I or K directions and for E , in the I or J directions. 

z 


If NOPE defines a two-dimensional surface, the appropriate E-field 

t sn 

components are set equal to the surface is in the X-Y plane, 

for example, E x on the two sides of the rectangle defined by the cell and E 

tan . 

on the other two sides are set equal to -E. . . ... 

^ incident 
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Figure 6. Determine E-Field Boundary Conditions 

For a lossy dielectric aircraft or an aircraft composed, in part, by 
a lossy dielectric, the boundary condition is given by 


a£ s . *s * 
e 3T + aE = ^ 


(e “ e o ) 



This boundary condition holds within the lossy dielectric volume. It may be 
considered a volume boundary condition. In the limit, o ■* «, it becomes the 
surface boundary condition 


E tan m _ £ tan 

scattered incident 

for perfectly conducting surfaces. 
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Solid blocks (3“dimensional bodies) and plates (2-dimensional bodies) 
are treated much the same as for the perfectly conducting counterparts, except: 

1) the volume boundary condition is used 

2) the NOPE array values of 1,2,3.12,13.23,123, and ^ are replaced by 
5,6,7,56,57,67,567 and 8 respectively 

3) the plates have thickness TX, TY and TZ in the x, y and z 
directions respectively, where the thickness is expressed as a 
fraction of the cell dimension in the appropriate direction 

*1) the volume boundary condition is slightly modified, as discussed 
discussed in Appendix A: Thin Plate Formalism, to more accurately 

model plates. 

The generalized flow of the Subroutine EBC is given in Figure 7. The 
details of the flow through the perfectly conducting portion of EBC are given 
in Figure 8. The lossy dielectric portion follows the same approach. 

Subroutine EBC is called from EADV. EINCX, EINCY, and EINCZ are called by EBC. 

Thin wires representing the lightning channels and interim wires (M=1 
and 12) are also input in EBC. 
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Figure 7. Generalized Flow Chart of Subroutine EBC 




Figure 8. Flow Chart of Subroutine EBC 
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SUBROUTINE EBC 


C 

COMMON/ EFIELD/EX( 29 ,29,29),EY(29,29,29) ,EZ(29,29,29) 
COMMON/HFIELD/HX(29,29,29) ,HY (29, 29, 29), HZ (29, 29, 29) 

COMMON/ GRID/ X( 28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(23) 

COMMON/ IS ITEM/ NOPE (29 ,29,29) 

COMMON/ EXTRAS/NX, NY, NZ ,NX1 ,NY1 ,NZ1 ,N,M,MQ,DT,XMU,EPSO,EPS,NPTLIM, 

1 NN ,NPTS,LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX.TY.TZ, AMP, ALPHA, BETA, IDLS 
C 

IF( IDLS.EQ.O) GO TO 1 
IF(M.NE.l) GO TO 2 
C 

UNEXPANDED LIGHTNING CHANNELS 

DO 3 J=13 , 28 
EY(8,J,15)=0.0 

3 CONTINUE 

DO A J-1,11 
EY(22,J,15)=0.0 

4 CONTINUE 

GO TO I 

2 CONTINUE 

IF(M.NE .2) GO TO 5 

EXPANDED LIGHTNING CHANNELS 

DO 6 J=9,28 
EY(5,J,13)=0.0 
6 CONTINUE 

GO TO 1 

5 CONTINUE 


I CONTINUE 


LOOP OVER REGION IN WHICH EBC MUST BE SET 
C 

DO 500 1=1, NX1 
DO 500 J-l.NYl 
DO 500 K=1 ,NZ1 
C 

IF(M.EQ.l) GO TO 12 

TX=0,008 

TY=0.02 

TZ=0.02 
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12 CONTINUE 


C 

IF(NOPE(I,J,K).EQ.O)GO TO 500 
C 

IF (N0PE(I,J,K) .NE.4) GO TO 60 
C 

C SET BOUNDARY CONDITIONS ON PERFECTLY CONDUCTING 3-D BODY 
C 

C SET EX 

C 

IF((J~1) .EQ.O.OR.(K-L).EQ.O) GO TO 18 

IF ( HOPE ( I , J , K- 1 ) . EQ. 4 .AND. NOPE ( I , J-l ,K~1 ) . EQ. 4 . AND. 

1 N0PE(I,J-1,K).EQ.4)G0 TO 10 

EX(I,J,K)=-EINCX(I,J,K) 

10 CONTINUE 

IF(NOPE(I,J+l,K) .EQ.4)G0 TO 15 
EX( I , J+l , K+l ) =~EINCX ( I , J+l , K+l ) 

EX( I , J+l , K)=-EINCX( I , J+l ,K) 

15 IF(NOPE(I,J,K+l).EQ.4)GO TO 20 
EX(I,J,K+1)=-EINCX(I,J,K+1) 

EX( I , J+l ,K+1)=-EINCX( I, J+l ,K+1) 

C 

SET EY 
18 CONTINUE 

IF((I~1) .EQ.O.OR.(K-l) .EQ.O) GO TO 38 
20 IF ( HOPE ( 1-1 , J ,K) .EQ.4.AND.N0PE(I-1,J,K-1) .EQ.4.AND. 

1 N0PE(I,J,K-1) .EQ.4)G0 TO 30 

EY(I,J,K)=-EINCY(I,J,K) 

30 CONTINUE 

IF(N0PE(I+1,J,K).EQ.4)G0 TO 35 
EY( 1+1 , J , K+l ) =-EINCY ( 1+1 , J , K+l ) 
EY(I+1,J,K)=-EINCY(I+1,J,K) 

35 IF(NOPE( I , J ,K+1) . EQ.4)G0 TO 40 
EY( I , J ,K+1 )=— EINCY( I , J , K+l ) 

EY( 1+1 , J ,K+1)=-EINCY ( 1+1 , J ,K+1 ) 

SET EZ 

38 CONTINUE 

IF( ( I— 1 ) .EQ.O. OR. (J-l) .EQ.O) GO TO 500 
40 IF(N0PE(I-1 ,J,K) .EQ.4.AND.H0PE(I-1 > J-1,K) .EQ.4.AND. 

1 NOPE (I, J-l ,K) .EQ.4)G0 TO 50 
EZ( I , J ,K)=— EINCZ( I, J ,K) 

50 CONTINUE 

IF(N0PE(I+1,J,K).EQ.4)G0 TO 55 
EZ(I+1,J,K)=-EINCZ(I+1,J,K) 

EZ( 1+1 , J+l ,K)=-EINCZ( 1+1 , J+l ,K) 

55 IF ( NOPE ( I , J+l ,K) .EQ. 4)G0 TO 500 
EZ( 1+1 , J+l ,K)=-EINCZ( 1+1 , J+l ,K) 

EZ( I, J+l ,K)=-EINCZ( I , J+l ,K) 

GO TO 500 

60 CONTINUE 
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IF(N0PE(I,J,K).NE.1.AND.N0PE(I,J,K) .NE. 12.AND.N0PE(I, J,K) .NE. 13. 

1 AND. HOPE ( I,J,K).NE. 123 .AND.NOPE( I , J ,K) .NE.2.AND. 

2 NOPE(I,J,K).NE.23.AND.NOPE(I,J,K).NE.3)GO TO 510 
C 

C SET BOUNDARY CONDITIONS ON PERFECTLY CONDUCTING 2-D BODY 
C 

IF(N0PE( I , J , K) .NE . 1 . AND. NOPE ( I , J , K) ,NE. 12 . AND. N0PE( I, J ,K) .NE . 13 . 
1 AND.NOPE( I, J ,K) .HE. 123) GO TO 70 

Y-Z PLANE 

EY( I , J , K+l ) =-E INCY ( I , J , K+l ) 

EY(I,J,K)=-EINCY(I,J,K) 

EZ( I , J+l , K) — E INCZ( I , J+l , K) 

EZ( I, J ,K)=-EINCZ( I, J ,K) 

IF( NOPE( I , J , K) .EQ. 12 .OR. NOPE ( I ,J,K).EQ.123) GO TO 70 
IF(N0PE( I, J ,K) .EQ. 13) GO TO 80 
GO TO 500 

70 IF(N0PE( I, J ,K) .NE. 2. AND.N0PE( I , J,K) .NE. 23. AND. 

1 NOPE( I, J,K) .NE. 12.AND.N0PE(I, J ,K) .NE.123) GO TO 80 

X-Z PLANE 

EX( I , J ,K+1 )=-EINCX( I , J , K+l ) 

EX( I, J ,K)=-EINCX( I , J ,K) 

EZ(I+1,J,K)=-EINCZ(I+1,J,K) 

EZ(I,J,K)*-EINCZ(I,J,K) 

IF(N0PE(I,J,K).EQ.23.0R.N0PE(I,J,K).EQ.123) GO TO 80 
GO TO 500 

80 CONTINUE 

X-Y PLANE 

EX( I , J+l ,K)=— EINCX( I , J+l ,K) 

EX(I,J,K)»-EINCX(I,J,K) 

EY( 1+1 , J ,K)=-EINCY( 1+1 , J,K) 

EY(I,J,K)=-EINCY(I,J,K) 

GO TO 500 

510 CONTINUE 

SET COMPUTATIONAL PARAMETERS 

8 RSIGMA-1. /SIGMA 

EPEFFX=(1.0-TX)*EPS0+(TX*EPS) 

EPEFFY=( 1 . 0-TY)*EPSOf (TY*EPS) 

EPEFFZ=( 1 ,0-TZ)*EPSO+(TZ*EPS) 

EXPON=EXP( (-SIGMA*DT) /EPS) 

EXPDX=EXP( (-SIGMA*TX*DT) /EPEFFX) 

EXPDY=EXP( (-SIGMA*TY*DT) /EPEFFY) 
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EXPDZ=EXP((-SIGMA*TZ*DT)/EPEFFZ) 

XEXP=1 . -EXPON 
XEXPX-l.-EXPDX 
XEXPY=1 .-EXPDY 
XEXPZ=1 .-EXPDZ 
C 

IF ( NOPE ( I , J ,K) .NE. 8)G0 TO 600 
C 

SET BOUNDARY CONDITIONS ON LOSSY DIELECTRIC 3-D BODY 
SET EX 

IF((J-1) .EQ.O.OR.(K-l) .EQ.O) GO TO 200 
EX(I,J,K)=(EX(I,J,K)*EXPON)+((XEXP)*(-EINCX(I,J,K) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCX(I, J,K) 

2 +(( (HZ( I, J,K)-HZ( I , J-l ,K) )*DYI( J) ) 

3 -( (HY( I, J ,K)-HY(I, J ,K-1 ) )*DZI(K) ) )*RSIGMA) ) 

IF(N0PE(I,J+1,K) .EQ.8)G0 TO 150 

EX(I, J+l ,K)=(EX(I , J+l ,K)*EXP0N)+( (XEXP)*(-EINCX(I, J+l ,K) 

1 -( (EPS-EPSO)*RSIGMA) *DEINCX( I , J+l ,K) 

2 +( ( (HZ(I , J+l ,K)-HZ( I , J,K) )*DYI( J+l) ) 

3 -( (HY( I , J+l ,K)-HY( I, J+l ,K-1) )*DZI(K)) )*RSIGtlA) ) 

EX( I , J+l ,K+1 ) =(EX( I , J+l ,K+1 )*EXP0N)+( (XEXP)*(-EINCX(I, J+l ,K+1 ) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCX( I , J+l ,K+1) 

2 +( ( (HZ(I , J+l ,K+1)-HZ(I , J,K+1) )*DYI( J+l) ) 

3 -( (11Y( I, J+l ,K+1)-HY( I, J+l ,K) )*DZI(K+1) ) )*RSIGMA) ) 

150 IF(N0PE(I,J,K+1) .EQ.8)G0 TO 200 

EX( I , J , K+l )=(EX( I , J , K+l ) *EXP0N )+( (XEXP) * ( -EINCX( I , J , K+l ) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCX( I, J ,K+1) 

2 +(((HZ(I,J,K+1)-IIZ(I > J-1,K+1))*DYI(J)) 

3 -( (HY( I, J,K+1)-HY( I , J,K) )*DZI(K+1) ) )*RSIGMA) ) 

EX( I , J+l ,K+1 ) =(EX( I , J+l , K+l )*EXP0N)+( (XEXP) *(-EINCX( I , J+l , K+l ) 

1 -( (EPS~EPSO)*RSIGMA)*DEINCX( I , J+l , K+l) 

2 +( ( (HZ( I , J+l ,K+1 )— HZ( I , J,K+1 ) )*DYI( J+l ) ) 

3 -( (HY(I, J+l ,K+1 )-HY( I, J+l ,K) )*DZI(K+1) ) )*RSIGMA) ) 

SET EY 
C 

200 CONTINUE 
C 

IF((I-1) .EQ.O. OR. (K-l) .EQ.O) GO TO 300 
EY(I,J,K)=(EY(I,J,K)*EXPON)+((XEXP)*(-EINCY(I,J,K) 

1 -( (EPS-EPSO) *RSIGMA)*DEINCY( I , J, K) 

2 +(((HX(I,J,K)-HX(I > J,K-1))*DZI(K)) 

3 -( ( HZ ( I , J , K)-HZ ( 1-1 , J , K) ) *DXI ( I ) ) ) *RSIGMA) ) 

C 

IF( NOPE( 1+1 , J , K) .EQ.8)GO TO 250 
C 

EY( 1+1 , J , K+l )= (EY( 1+1 , J , K+l ) *EXPON )+( (XEXP) * (-EINCY( 1+1 , J , K+l ) 
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1 -((EPS-EPS0)*RSIGHA)*DEINCY(I+1,J,K+1) 

2 +(((HX(I+1,J,K+1)-HX(I+1,J,K))*DZI(K+1)) 

3 -((HZ(I+1,J,K+1)-HZ(I,J,K+1))*DXI(I+1)))*RSIGMA)) 

C 

EY( 1+1 , J,K)=(EY( 1+1 , J ,K)*EXPON)+( (XEXP)*(-EINCY(I+1 , J,K) 

1 -( ( EPS-EPSO ) *RSIGMA) * DE INCY ( 1+1 , J , K) 

2 +(((HX(I+1,J,K)-HX(I+1,J,K-1))*DZI(K)) 

3 -( (HZ( 1+1 , J ,K)-HZ( I, J,K) )*DXI( 1+1) ) )*RSIGMA) ) 

C 

250 IF( NOPE ( I , J , K+l ) . EQ . 8 ) CO TO 300 
C 

EY(I,J,K+1)=(EY(I,J,K+1)*EXP0N)+((XEXP)*(-EINCY(I,J,K+1) 

1 -((EPS-EPS0)*RSIGMA)*DEINCY(I,J,K+1) 

2 +(((HX(I,J,K+1)-HX(I,J,K))*DZI(K+1)) 

3 — ( (11Z( I, J ,K+1 )-HZ( 1-1 , J ,K+1) )*DXI( I) ) )*RSIQ'1A) ) 

C 

EY(I+1,J,K+1) =(EY(I+1,J,K+1)*EXP0N)+((XEXP)*(-EINCY(I+1,J,K+1) 

1 -(( EPS-EPSO )*RSIGMA)*DEINCY( 1+1 , J ,K+1) 

2 +(((HX(I+1,J,K+1)-HX(I+1,J,K))*DZI(K+1)) 

3 -( (HZ(I+1 , J,K+1)-HZ( I, J,K+1) )*DXI(I+1) ))*RSIGMA) ) 

SET EZ 

300 CONTINUE 
C 

IF((I-I) .EQ.O.OR. (J-l) .EQ.O) GO TO 500 
EZ(I,J,K)*(EZ( I , J,K)*EXPON)+( (XEXP)*(-EINCZ( I, J ,K) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCZ( I, J, K) 

2 +(((HY(I,J,K)-HY(I-1,J,K))*DXI(I)) 

3 -((HX(I,J,K)-HX(I,J-1,K))*DYI(J)))*RSIGMA)) 

C 

IF(N0PE(I+I,J > K).EQ.8)G0 TO 350 
C 

EZ(I+1,J,K)=(EZ(I+1,J,K)*EXP0N)+((XEXP)*(-EINCZ(I+1,J,K) 

1 -( (EPS-EPSO) *RSIG11A) *DEINCZ( I+I , J,K) 

2 +(((HY(I+I,J,K)-HY(I,J,K))*DXI(I+1)) 

3 -( (HX( I+I , J ,K)-HX( 1+1 , J-l ,K))*DYI(J) ) )*RSIGMA) ) 

C 

EZ(I+1 , J+l ,K) =(EZ( I+I , J+I ,K)*EXPON)+( (XEXP)*(-EINCZ( 1+1 , J+l ,K) 

1 -((EPS-EPS0)*RSIGMA)*DEINCZ(I+1,J+1,K) 

2 +( ( (HY( 1+1 , J+l , K)-IIY( I , J+l , K) ) *DXI ( 1+1 ) ) 

3 -((UX( 1+1 , J+l ,K)-HX( 1+1 , J ,K) )*DYI( J+l) ))*RSIGHA) ) 

C 

350 IF(N0PE(I,J+1,K).EQ.8)G0 TO 500 
C 

EZ(I+1 , J+l ,K) =(EZ( 1+1 , J+l ,K)*EXPON)+( (XEXP)*(-EINCZ(I+1 , J+l ,K) 

1 -( (EPS-EPSO) *RSIGMA) *DEINCZ( 1+1 , J+l , K) 

2 +( ( (HY( 1+1 , J+l , K)-HY( I , J+l , K) ) *DXI ( 1+1 ) ) 

3 - ( ( HX( 1+1 , J+l , K)-HX( 1+1 , J , K) ) *DYI ( J+l ) ) ) *RS IGMA) ) 

C 

EZ( I , J+l ,K)=(EZ( I , J+l ,K)*EXPON)+( (XEXP)*(-EINCZ( I , J+l ,K) 
l-( ( EPS-EPSO )*RSIGMA)*DEINCZ( I, J+l ,K) 

2 +( ( (HY(I , J+l ,K)-HY(I-1 , J+l ,K) )*DX1( I) ) 

3 -((HX(I,J+1,K)-HX(I,J,K))*DYI(J+1)))*RSI01A)) 

C 


77 


non no 


GO TO 500 


C 

600 CONTINUE 
C 

SET BOUNDARY CONDITIONS ON LOSSY DIELECTRIC 2-D BODY 

IF((I-I) .EQ.O.OR.(J-l) .EQ.O.OR.(K-l) .EQ.O) GO TO 500 
IF ( NOPE ( I, J,K) .NE.5.AND.HOPE( I, J,K) .NE.56.AND.NOPE(I, J,K) . 

1 NE . 57 . AND.NOPE(I , J ,K) .NE. 567) GO TO 700 

Y-Z PLANE 

EY(I,J,K+1)=(EY(I,J,K+1)*EXPDX)+((XEXPX)*(-EINCY(I,J,K+1) 

1 -( (EPS-EPSO )*RSIGMA)*DEINCY( I ,J ,K+1) 

2 +( ( (HX( I , J , K+l )-HX( I , J ,K) )*DZI( K+l) ) 

3 -( (HZ( I , J ,K+1 )-HZ( I-I , J ,K+1 ) )*DXI( I) ) )*RSIGMA*TX) ) 

C 

EY( I , J ,K)=(EY( I , J , K)*EXPDX)+( (XEXPX)*(-EINCY(I , J ,K) 

1 -((EPS-EPSO)*RSIGMA)*DEINCY(I,J,K) 

2 +( ( (HX( I , J ,K)-HX( I , J , K-l ) )*DZI(K) ) 

3 -( (HZ( I , J,K)-HZ( 1-1 , J ,K) )*DXI( I) ) )*RSIGMA*TX) ) 

C 

EZ(I,J+1,K)=(EZ(I,J+I,K)*EXPDX)+((XEXPX)*(-EINCZ(I,J+1,K) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCZ( I, J+l ,K) 

2 +( ( (HY(I, J+l ,K)-HY(I-1 , J+l ,K) )*DXI(I) ) 

3 -( (HX( I , J+l ,K)-tIX( I , J ,K) )*DYI( J+l) ) )*RSIGMA*TX) ) 

C 

EZ( I , J , K)-(EZ( I , J , K)*EXPDX)+( (XEXPX) *(-EINCZ( I , J ,K) 

1 -( (EPS _ EPSO)*RSIGMA)*DEINCZ( I, J,K) 

2 +( ( (HY( I , J ,K)-HY(I-1 , J,K) )*DXI( I) ) 

3 -( ( HX( I , J ,K)-HX( I , J— 1 ,K) )*DYI( J) ) )*RSIGMA*TX) ) 

C 

IF(NOPE(I,J,K).EQ.56) GO TO 700 
IF(NOPE(I,J,K) .EQ.57) GO TO 800 
IF(N0PE(I, J,K) .EQ.567) GO TO 700 
GO TO 500 
C 

700 IF(NOPE(I,J,K) . NE. 6 . AND. NOPE ( I , J ,K) .NE. 67 . AND. 

1 N0PE(I,J,K) .NE.56.AND.NOPE(I,J,K) .NE.567) GO TO 800 
C 

C X-Z PLANE 

C 

EX(I,J,K)=(EX(I,J,K)*EXPDY)+((XEXPY)*(-EINCX(I,J > K) 

1 -( (EPS-EPSO)*RSIGMA) *DEINCX( I , J ,K) 

2 +( ( (HZ( I , J ,K)-11Z( I , J-l ,K) )*DYI( J) ) 

3 -( (HY( I , J ,K)-HY( I , J ,K-1 ) )*DZI(K) ) )*RSIGMA*TY) ) 

C 

EX(I,J,K+1)=(EX(I,J,K+1)*EXPDY)+((XEXPY)*(-EINCX(I,J,K+1) 

1 -( (EPS-EPSO) *RSIGilA)*DEINCX( I , J ,K+1) 

2 +( ( (HZ( I , J ,K+1 )-HZ( I , J-l ,K+1 ) )*DYI( J) ) 

3 -((liY(I,J,K+l)-HY(I,J,K))*DZI(K+l)))*RSIGMA*TY)) 
C 

EZ( 1+1 , J ,K)=( EZ( 1+1 , J ,K)*EXPDY)+( (XEXPY) *(-EINCZ(I+l , J ,K) 

1 -( (EPS-EPSO) *RSIGMA)*DEINCZ( 1+1 , J,K) 

2 +( ( (HY(I+1 , J,K)-UY(I, J , K) )*DXI(I+1) ) 
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C 3 -( (HX( 1+1 , J ,K)-HX( 1+1 , J-l ,K) )*DYI( J) ) )*RSIGMA*TY) ) 

EZ(I,J,K)=(EZ(I,J,K)*EXPDY)+((XEXPY)*(-EINCZ(I,J,K) 

1 -((EPS-EPSO)*RSIGMA)*DEINCZ(I,J,K) 

2 +(((HY(I, J,K)-HY(I-1 , J,K))*DXI(I)) 

3 -( (HX( I, J ,K)-HX( I, J-l ,K) )*DYI( J) ) )*RSIGI1A*TY) ) 

C 

IF(NOPE(I, J,K) .EQ.67 . OR.NOPE(I, J,K) .EQ.5G7) GO TO 800 
GO TO 500 

800 CONTINUE 

X-Y PLANE 

EX( I , J ,K)=(EX( I , J ,K)*EXPDZ)+( (XEXPZ) *(-EINCX( I , J ,K) 

1 -((EPS-EPSO)*RSIGMA)*DEINCX(I,J,K) 

2 +(((HZ(I,J,K)-IIZ(I,J-1,K))*DYI(J)) 

3 -( (HY( I , J ,K)-HY( I,J,K-1))*DZI(K) ) )*RSIGMA*TZ) ) 

EX( I , J+I , K) =( EX( I , J+l , K) *EXPDZ )+ ( ( XEXPZ )*(-EINCX( I, J+l , K) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCX( I , J+l ,K) 

2 +( ( (HZ( I, J+I ,K)-I1Z( I , J ,K) )*DYI( J+l ) ) 

3 -((UY(I,J+l,K)-tiY(I,J+l,K-l))*DZI(K)))*RSIGMA*TZ)) 
EY(I,J,K)=(EY(I,J,K)*EXPDZ)+((XEXPZ)*(-EINCY(I,J,K) 

1 -( (EPS-EPS0)*RSIGMA)*DEINCY( I , J,K) 

2 +( ( ( liX ( I , J ,K)— HX( I , J ,K-1 ) )*DZI(K) ) 

3 -( (HZ( I, J ,K)-HZ( I-I , J,K) )*DXI( I) ) )*RSICMA*TZ) ) 

EY( 1+1 , J,K)=(EY(I+1 , J ,K)*EXPDZ)+( (XEXPZ) *(-EINCY( I+I , J,K) 

1 -( (EPS-EPSO)*RSIGMA)*DEINCY( 1+1 , J ,K) 

2 +( ( (HX( 1+1 , J,K)-HX( 1+1 , J , K-l ) )*DZI(K) ) 

3 -( (UZ(I+I , J ,K)-HZ( I , J,K) )*LXI(I+1) ) )*RSIGMA*TZ) ) 

500 CONTINUE 

501 CONTINUE 

NWIRE=1 

IF(NWIRE.EQ.O) GO TO 550 
IF(M.NE.2) GO TO 550 

B1 TO B2 

1=4 
K.= 13 

DO 900 J=10, 14 
EY(I,J,K)=-EINCY(I,J,K) 

900 CONTINUE 
C 

EX(4 , 10, 13)=-EINCX(4 , 10,13) 

EZ(5, 10, 13)=-EINCZ(5, 10,13) 

EZ( 5 , 10 , 14)=-EINCZ(5 , 10,14) 

EX(5 , 10, 15)=-EINCX(5 , 10 , 15) 

EZ(6 , 10, 15)=-EINCZ(6 , 10,15) 

EZ(6, 10 , 16)=-EINCZ(6 , 10, 16) 
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B2 THROUGH B4 

EX( 6 , 10 , 1 7)=-EINCX( 6,10,17) 
EZ(7 , 10, 17)=-EINCZ(7 ,10,17) 


J=10 
K=18 

DO 901 1=7,10 
EX( I , J , K)=-EINCX( I , J , K) 

901 CONTINUE 

AFT OF B4 THROUGH B6 

EZ(11 , 10, 18)=-EINCZ( 11,10,18) 

J=10 
K=19 

DO 902 1=11,16 
EX(I,J,K)=-EINCX(I,J,K) 

902 CONTINUE 

AFT OF B6 TO B8 

1=17 
K=19 

DO 903 J=10, 13 
EY(I,J,K)=-EINCY(I,J,K) 

903 CONTINUE 

EX( 17 , 14 , 19 )=-EINCX( 17,14,19) 
EY(18, 13, 19)=-EINCY(18 , 13,19) 
EX(18, 13, 19)=-EINCX( 18 , 13,19) 
LY(19,13, 19)=-EINCY(19 , 13 , 19) 
EX( 19,14, 19)=-EINCX( 19,14,19) 

B8 ON TO END 

EX(20, 14 , 19)=-EINCX(20, 14,19) 
EZ(21,14,18)=-EINCZ(21,14,18) 
EZ(21 , 14, 17)=-EINCZ(21 ,14,17) 

J=14 
K=17 

DO 904 1=21,25 
EX( I , J , K)=— EINCX( I , J , K) 

904 CONTINUE 
C 

EY(26, 13 , 17)=-EINCY(26, 13 , 17) 
EY(26, 12, 17)=-EINCY(26 , 12,17) 
EY(26, 11 , 17)=-EINCY(26 , 11 ,17) 
C 

EX(26 , 11 , 17)=-EINCX(26 ,11,17) 
C 

EY(27, 11 , 17)=-EINCY(2 7 , 11 , 17) 
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EY(27 , 12 , 17)=-EINCY(27, 12,17) 
C 

EX(27 , 13, 17)=-EINCX(27 ,13,17) 
C 

550 CONTINUE 
C 

RETURN 

C 

END 
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3.5.5 


Subroutine ERAD 


Subroutine ERAD evaluates the expression 

E (R,D,EL1 ,EL2,EL3,L1 ,L2,L3,TH) * R/D* (1) 

(ELI * (L2 - L3) * (L2*L3 + TH * (L2 + L3) + TH**2) + 

EL2 * (L3 “ LI) * (L3*L1 + TH * (L3 + LI ) + TH**2) + 

EL3 * (LI - L2) * (LI *L2 + TH * (LI + L2) + TH**2)). 

It substitutes the results for the EADV-determined values of the field 
components which are on the surface of the problem cell space (shown in 
Figure 9) according to the prescription: 

EX(1, 1 ,K) = E (RXY(I,K,1) ,D,EX(I,2,K) , EXYD(I,K,2), 

EXYD(I,K,3) ,L1 ,L2,L3, THXY, ( I , K, 1 ) ) (2a) 


EZ(I,NY,K) = E ( RZY ( I , K, 2) , D, EZ (I ,NY1 ,K) , EZYU(I,K,2), 

EZYU(I,K,3),L1 ,L2,L3,THZY(I,K,2)) (21) 

This substitution imposes a far field behavior on the outermost 
region of the problem cell space. Equation (1) is derived directly from the 
assumption that the outermost electric field components of Figure 9 and the 
corresponding electric fields that are one cell in from the outer surface, 
behave as far fields,* that is, 

E - f(e,<t>)g(t-r/c)/r or equivalently (3) 

Er - f(e,<t>)g(t-r/c) 

*As a result, fields on the outermost edges are not written over, but this does 
not effect the interior results. 
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Figure 9. Field 



at the Outer Boundary of the Cell Space 


Thus, taking E^ +1/,2 (I,J,K) as an example, 


R (1,J,K) e" + 1/2 (1,J,K) = e ( J, K, [ (n+1 /2) At - R ( 1 ,J,K)/c]) (4) 

ja y y yx 

R (2,J,K)e" +1/2 (2,J,K) s e ( J, K, [ (n+1 /2) At - R(2 , J, K) /c] ) (5) 

yx y y yx 

where R is the distance to the E component, 
yx y 

What is desired is E^ +1 ^ 2 (I , J, K) expressed in terms of E^ +1/ ^ 2 *1 

(2,J,K), i.e., the outer boundary field component expressed in terms of an 
interior corresponding field component at various time offsets(2,^) . Equation 
(1) is, therefore, rewritten as, 


R (1 ,J,K) E n+1/2 (1 ,J,K) = e 
yx y y 


R (2 , J , K) R ( 1 , J , K) - 

yx yx * 


V 2 > J » K) j^ 


1/2) At + At 


= e 


(j,K, £(n- 


R (2.J.K) 

i /2) At — ^ + e (j 

c yx 


,K,)AtJ ^ 


( 6 ) 


where 


R ( 1 , J, K) - R (2 , J, K) 

0 5 i _ yx 

yx cAt 


(7) 


and Equation (2) is rewritten with a time shift as, 


R(2 , J, K) Ey~ 1/2 U (2,J,K) 


(j,K, £(n- 


R (2, J,K) 
1/2) At - ~ 


- ijAtJ) 


(8) 
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It can be seen that Equations (6) and (8) are identical except for t i 

in place of 0 in Equation (8). Either Equation (6) or (8) can be expanded in 

terms of the time 14 t or 6 At. The first three terms of Maclaurin's series 

i yx 

are used so that 

R ( 1 , J , K) E n+1/2 ( 1 , J , K) * A + B (0 At) + C (0 At) 2 (9) 

yx y yx yx 

R yx (2,J,K) Ey _1 /2_ ^ i (2 > J,K) = A + B (H^At) + C (X^At) 2 . (10) 

If three different values of are selected, then, from 

Equation (10), we obtain three equations for the three unknown A, B and C in 

terms of the three values of E n_1 /2 4 i, the three values of % and R (2,J,K), 

y l y a 

all of which are know if JL i -1 . Therefore, A, B and C can be determined and 
are 

A - [Ey~ 1/2 ~ !li (2,J,K)jl 2 t 3 0l 2 - l 3 ) (ID 

+ E y -1 /2-2 -i(2 , J , K) (& 3 - ) 

+ Ey~ 1/2 ~ Zi (2,J,K)Z 1 Z 2 U 1 - 1 2 )] R yx (2,J,K)/D 

B * [E y “ 1/2-itl (2,J,K) U 2 2 - l 3 2 ) (12) 

+ E y ~ 1/2-!l3 (2,J,K) (J^ 2 - t 2 2 )] R yx (2,J,K)/DAt 

C = [E n “ 1/2 ~ S ’ 1 (2,J,K) (i - l ) (13) 

y $ 

+ E y ~ 1/2_!! '2( 2> j f K) U 3 - t 1 ) 

+ E n ~ U2 ~ l H2,J,K) (JL - £.)] R (2 , J , K) /DAt 2 
y * ^ yx 
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where D is the determinant 


d - i 2 i 3 a 2 - * 3 ) ♦ + a 1 i 2 u 1 - % 2 ) 

Knowing A, B and C, Equation (9) immediately yields 


( 14 ) 


, A + B(0 At) + C(0 At)‘ 

e" +1/2 (i,j,k) = yx yx 


R ( 1 , J , K) 
yx 


(15) 


This is the desired result of E^ +1 /2 (I , J, K) expressed in terms of E^ 1/2 *' i 

and other known quantities, namely R (I,J,K), R (2,J,K) and 8., 8. and 8_ 

yx yx i d. j 


The last equation, Equation (15), is identical to the expression. 

Equation (1), that ERAD evaluates. 

For example, the terms EX(I,2,K), EXYD(I,K,2), and EXYD(I,K,3) in 

Equation (2a) are the backstored values of EX at J = 2 that are used to 

calculate EX(I,1,K). They are found in Subroutine EADV. The quantities D, LI, 

L2, L3, THXY(I,K, 1) correspond to D, 8., 8 , and 0 (I,K). D is calculated in 

id x y 

ERAD as are 8 1 , 8 2 and 8^ and the corresponding stored values of E one cell in 

from the outer boundary. The 8^'s are spread enough in time to correspond to 

the propagation time across the expanded cell boundary. Reference 2 details 

what requirements on the 8., and, hence, on N2 and N3 in EADV must be met. Only 

when constant mesh size is used can the 8 1 = -1 , 8 2 =0, 8^=1. In this case, the 

expression in Equation (1) goes over exactly to the expression in Section 2.3. 

The value of R (I,K,1) (which represents R (I,K,2)/R (I,K,1) in this 

xy xy xy 

subroutine) is determined in Subroutine SETUP, as is the value of THXY (I,K, 1). 
Since they are geometry dependent quantities, they need only be found once, 
hence, their presence in SETUP. 


Flow for Subroutine ERAD is as shown in Figure 10. Local variables 
are described in Table 5. Subroutine ERAD is called by EADV. 
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Table 5. Local Variables - Subroutine ERAD 


Variable 

NAME 


D 

LI 

L2 

L3 



Determinant 

Time shifts for backstored data in 
number of cycles 
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Start 



Figure 10. Flow Chart for Subroutine ERAD 
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SUBROUTINE ERAD 


C0MM0N/EFIELD/EX(29,29,29) ,EY(29,29,29) ,EZ(29,29,29) 
C0MM0N/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX, NY, NZ ,NX1 ,NY1 , NZ1 ,N,M,MQ,DT, XMU.EPSO, EPS, NPTLIM, 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, 

2 BETA,IDLS 

COMMON/RAD/RYX(28,29,2) ,RZX(29, 28, 2) ,RXY(28 , 29,2) ,RZY(29,28,2) , 

1 RXZ(28,29,2) ,RYZ(29,28,2) ,THYX(28,29,2) ,TUZX(29,28,2) , 

2 THXY(28,29,2) ,THZY(29,28,2) ,THXZ(28,29,2) ,THYZ(29,28,2) 
COMMON/EBS/EYXD(28,29,3),EYXU(28,29,3),EZXD(29,28,3),EZXU(29,28,3) 

1 ,EXYD(28,29,3) ,EXYU(28 , 29 , 3) ,EZYD(29,28,3) ,EZYU(29,28,3) , 

2 EXZD(28,29,3) ,EXZU(28,29,3) ,EYZD(29 ,28 , 3) ,EYZU(29 ,28, 3) , 

3 N1,N2,N3 

E(R,D,EL1,EL2,EL3,L1,L2,L3,TH) = R/D*( 

1 EL1*(L2-L3)*(L2*L3+TH*(L2+L3)+TH**2)+ 

2 EL2* (L3-L1 ) *(L3*L1+TH* ( L3+L1 )+TH**2 )+ 

3 EL3*(L1-L2)*(L1*L2+TH*(L1+L2)+TH**2)) 


LI— 1 

L2=N-N2-1 

L3=N-N3-1 

D=L2*L3*(L2-L3)+L1*L3*(L3-L1)+L1*L2*(L1-L2) 

IF(D.EQ.O.) D=l.E-9 

RADIATE EXY 

DO 11 1=1, NX1 
DO II K-2.NZ1 

EX(I, 1 ,K)=E(RXY( I,K, 1) ,D,EX(I,2,K) ,EXYD(I,K,2) ,EXYD(I,K,3) , 

1 LI ,L2 ,L3 ,THXY( I,K, 1) ) 

EX( I , NY, K) =E(RXY( I,K, 2) ,D,EX( I, NY1 ,K) ,EXYU(I,K,2) ,EXYU(I,K,3) , 
1 LI ,L2 ,L3 ,THXY( I,K, 2) ) 

II CONTINUE 

RADIATE EXZ 

DO 111 1=1, NX1 
DO 111 J=2 ,NY1 

EX(I, J , 1)=E(RXZ( I , J , 1) ,D,EX(I , J , 2) ,EXZD(I , J,2) ,EXZD(I, J ,3) , 

1 LI ,L2 ,L3 ,THXZ( I, J, 1)) 

EX(I, J,NZ) =E(RXZ(I, J,2) ,D,EX(I, J,NZ1) ,EXZU(I,J,2) ,EXZU(I, J,3) , 
1 LI ,L2,L3,THXZ(I,J,2)) 

111 CONTINUE 

RADIATE EYX 

DO 22 J=1 ,NY1 
DO 22 K=2 ,NZ1 

EY( 1 , J,K)=E(RYX( J ,K, 1) ,D,EY(2 , J , K) ,EYXD( J ,K, 2) ,EYXD(J,K,3) , 

1 LI ,L2 ,L3 ,THYX( J,K, 1) ) 
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EY(NX,J,K) =E(RYX(J,K,2),D,EY(NX1,J,K),EYXU(J,K,2),EYXU(J,K,3) , 
1 L1,L2,L3,THYX(J,K,2)) 

22 CONTINUE 
C 

C RADIATE EYZ 
C 

DO 222 1=2, NX1 
DO 222 J=1 ,NY1 

EY( I, J , 1)=E(RYZ( I, J ,1 ) ,D,EY(I,J,2) ,EYZD(I , J,2) ,EYZD(I, J ,3) , 

1 L I , L2 , L3 , TH YZ( I , J , 1 ) ) 

EY(I, J,NZ) =E(RYZ( I, J , 2) ,D,EY(I,J,NZ1) ,EYZU(I, J ,2) ,EYZU(I , J , 3) , 
1 LI ,L2 ,L3 ,THYZ( I, J,2) ) 

222 CONTINUE 
C 

C RADIATE EZX 
C 

DO 33 J=2,NY1 
DO 33 K=1,NZI 

EZ( 1 , J ,K)=E(RZX( J ,K, 1),D,EZ(2,J,K) ,EZXD( J ,K,2) ,EZXD( J ,K, 3) , 

1 LI ,L2 ,L3,THZX( J,K, 1)) 

EZ(NX,J,K) =E(RZX(J,K,2),D,EZ(NX1,J,K),EZXU(J,K,2) ,EZXU(J,K,3) , 
1 L1,L2,L3,THZX(J,K,2)) 

33 CONTINUE 
C 

C RADIATE EZY 
C 

DO 333 1=2, NX1 
DO 333 K=1 ,NZ1 

EZ( I , 1 ,K)=E(RZY( I ,K, 1) ,D,EZ(I,2,K) ,EZYD( I , K, 2) ,EZYD(I , K, 3) , 

1 L1,L2,L3,THZY(I,K,1)) 

EZ( I ,NY,K) =E(RZY(I,K,2) ,D,EZ( I ,NY1 ,K) ,EZYU(I,K,2) ,EZYU(I,K,3), 
1 L1,L2,L3,THZY(I,K,2)) 

333 CONTINUE 
C 

RETURN 

END 
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3.5.6 


Subroutine HADV 


In subroutine HADV, the equation 


V x E = 



( 16 ) 


is applied to the scattered fields and is used to fine H from E and H at 
earlier times. This "time stepping" solution is found for each component of H, 
i.e., Equation (16) is written explicitly for each component as 


3H S 3E S 3E® 

x z y 

-y = - 

3t 3y 3z 


3H S 3E S 3E S 

y x z 

3t 3z 3x 


3H 3E 3E 

Z X X 

~y = - 

3t 3x 3y 

These equations are cast in their finite-difference form and solved for H at 
the latest time. Thus 


n n-1 /e"" 1/2 (I,J+1,K) - e" 1/2 (I,J,K) 

H x ( I* J 'K) " d.J.K) - (At/y) (-* (J+1) ' . Z Y (J) 


(17) 


eJM /2 (I,j,k+1) _ E n-1/2 (I>J>K) 

Z (K+1) - Z (K) 
o o 


n n-1 / C 1/2(I ' J ' K+1) ~ e"" 1/2 (I,J,K) 

Hj(I,J,K) = H" ( I, J, K) - (At/y) ( -5 ——JL^ 


(18) 


E n_1 /2 (I+1 ,J,K) - E n_1 /2 (I,J,K)'' 
z z 

x (i+i) - xji) 

0 o 
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( 19 ) 



It is assumed that y = y Q for all regions of the problem space. 

HADV is called every program cycle by DRIVER to advance the H-fields, 
just as Subroutine EADV is called every program cycle to advance the E-fields. 
It uses the prior values of the fields to calculate the new values. Initially, 
these prior values are zero. The first nonzero values are those values for E. 
These values are accessed by a call to EBC in the Subroutine EADV for advancing 
the E-fields. Thus, on the next program cycle these values will be available 
and nonzero values of the fields off the boundary will start appearing when 
Equations (17), 08), and (19) are advanced again. Since EBC is always called 
after the fields are advanced, the proper boundary conditions are maintained. 

HADV also defines a loop of H fields that drives the lightning channel 

for M«1 . 


The only variable used locally in HADV is DTMU, which is equal to 
At/y. Flow for the subroutine is shown in Figure 11. 
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SUBROUTINE UADV 


C 

COMMON/EFIELD/EX(29 ,29,29) ,EY(29,29,29) ,EZ(29,29,29) 
COMMON/llFIELD/HX(29,29,29) ,HY (29, 29, 29), HZ (29, 29, 29) 

COMMON/ GRID/X( 28) ,Y(28) ,Z(28) ,XO(29) ,Y0(29) ,ZO(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NXI , NY 1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS .NPTLIM , 
1 NN, NPTS.LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
C 

DTXMU = DT/ XMU 

ADVANCE HX 

DO 1 I = 2, NXI 
DO 1 J = l.NYl 
DO 1 K = 1,NZ1 

HX(I,J,K) = HX( I , J , K) - DTXMU* ( (EZ( I , J+l , K) - EZ(I, J,K) )*DYOI(J) 

1 - (EY( I , J ,K+I) - EY( I , J ,K) )*DZOI(K) ) 

1 CONTINUE 

ADVANCE HY 

DO 2 I = 1,NX1 
DO 2 J = 2 ,NY1 
DO 2 K = 1 ,NZ1 

HY(I, J,K) = HY( I , J ,K) - DTXMU* ( (EX(I, J , K.+1) - EX(I, J.LC) )*DZOI(K) 

I - (EZ(I+1,J,K) - EZ( I, J ,K) )*DXOI(I) ) 

2 CONTINUE 

ADVANCE HZ 

DO 3 I = 1.NX1 
DO 3 J = l.NYl 
DO 3 K = 2.NZ1 

HZ(I,J,K) = HZ(I.J.K) - DTXMU*( (EY( 1+1 , J ,K) - EY(I, J,K) )*DXOI(I) 

1 - (EX( I , J+l ,K) - EX( I, J.K) ) *DY0I( J) ) 

3 CONTINUE 


IF(IDLS.EQ.O) GO TO 1 
C 

IF(M.NE.l) GO TO 1 
C 

C M=1 LOOP , H-LOOP SOURCE 
C 

C M=2 LOOP NOT TREATED-SOURCE MUST BE OUTSIDE THE PROBLEM SPACE 
C AND THE WIRES ARE ALSO OUTSIDE 
C 

EXPAB=EXP ( -ALPHA* T ) -EXP ( -BET A* T ) 

C 

HX(8 , 21 , 15 >=-(33333 . )*EXPAB 
HZ(8, 21, 15)=( 16667. )*EXPAB 
HX(8 ,21 , 14 >=(33333. )*EXPAB 
HZ (7, 21, 15) =-(16667 . >*EXPAB 
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1 CONTINUE 


RETURN 

END 
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3 . 5.7 


FIELDS 


There are nine functions* used for calculating the six incident E- and 
H-field components: EINCX, EINCY, EINCZ, EINCXP, EINCYP, EINCZP, HINCXP, 

HINCYP, and HINCZP. The first letter in the function, E or H, tells what 
field; the next three letters, INC, indicate that these are incident field 
components; the next letter, X, Y or Z, states what component is evaluated at a 
point in space located at XP, YP, ZP and if P is not present the field 
component is evaluated at the field component location in the cell space 
labeled by I,J,K. EINCX, EINCY, and EINCZ are called from Subroutine EBC, The 
others are called by Subroutine DATASAV. 

Two different sources for the incident fields are considered: 

1) direct strikes; 

2) indirect illumination 

Direct strikes are modeled by H-fields impressed around a perfectly 
conducting wire leading to the aircraft that represents the lightning channel. 
If a lightning channel resistance is desired, then use the method outlined in 
Appendix D, in the discussion on excitation sources. The specified H-fields 
determine the lightning channel's current time history. Indirect illumination 
is modeled by the fields directly illuminating the aircraft, such as EINCX and 
HINCZ broadside illumination from above. 


* While these functions are not subroutines per se, they are similar enough in 
nature to warrant inclusion in this section, as opposed to having a separate 
short section of their own. 
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FUNCTION EINCX(I, J,K) 

C 

COHMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 
COMMON/GRID/X(28) ,Y(28) ,Z(28) ,XO(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZOI(28) 
COMMON/EXTRAS/NX , N Y , NZ , NXI , N Y1 , NZ l , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 

1 NN, NPTS, LMAX, SIGMA, C, T, PI, EXPFAC, I P,TX,TY,TZ, AMP, ALPHA, BETA, IDLS 

THREAT 

EINCX =0.0 
IF(IDLS.EQ.l) GO TO 1 
TAU = T - (UY0(23) - Y0(J))/C 
IF (TAU. LE. 0.0) GO TO 1 

E INCX=AMP* ( EXP ( -ALPHA*TAU)-EXP (-BETA*TAU) ) 

1 CONTINUE 
RETURN 
END 

FUNCTION DEINCX( I , J ,K) 

COtlMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZO(29) 

COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,ZO(29) , 

1 DX(29) , DY ( 2 9 ) ,DZ(29) ,DX0(28) ,DYO(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX, NY, NZ , NXI, NY1 ,NZ1 ,N,M,MQ ,DT, XMU, EPSO, EPS, NPTLIM, 

1 NN,NPTS,LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
DEINCX=0 . 0 

IF(IDLS.EQ.l) GO TO I 
TAU=T-(UY0(23)-Y0(J))/C 
IF(TAU.LE.O) GO TO 1 

DEINCX=AMP*(BETA*(EXP(-BETA*TAU) )-ALPHA* ( EXP(-ALPHA*TAU) ) ) 

1 CONTINUE 
RETURN 
END 

FUNCTION EINCY(I,J,K) 

COMMON/UGRID/UX( 28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 

COMMON/ GRID/X(28) ,Y(28) , Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX , NY , NZ , NXI , N Y1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 

THREAT 

EINCY = 0.0 
RETURN 
END 

FUNCTION DEINCY(I,J,K) 

DEINCY=0.0 
RETURN 
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END 


C 

FUNCTION EINCZ(I,J,K) 

C 

COMMON/ UGRID/UX( 28) ,UY(28) ,UZ(28) ,UX0(29) ,UYO(29) ,UZ0(29) 
C0MM0N/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 
I NN ,NPTS , LMAX, SIGMA,C ,T , PI , EXPFAC , IP ,TX,TY ,TZ , AMP, ALPHA, BETA,IDLS 

THREAT 


EINCZ = 0.0 

RETURN 

END 

FUNCTION DE INCZ( I , J ,K) 

COMMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 
COMMON/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) , YO(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX, NY, NZ ,NX1 ,NY1 ,NZ1 ,N,M,MQ,DT, XMU, EPSO, EPS, NPTLIM, 
I NN.NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
DEINCZ=0.0 
RETURN 
END 

FUNCTION HINCX(I, J,K) 

C0MM0N/UGRID/UX( 28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 
C0MM0N/GRID/X(28),Y(28),Z(28),X0(29),Y0(29),Z0(29), 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY1 , NZ 1 , N , M , MQ , DT , XMU , EPSO , EPS , NPTLIM , 
1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 

THREAT 

HINCX=0.0 

RETURN 

END 

FUNCTION HINCY ( I , J , K) 

C OMMON / UGRID/ UX( 2 8 ) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 

COMMON/ GRID/ X( 28 ) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY1 , NZ 1 , N , M ,MQ , DT , XMU , EPSO , EPS , NPTLIM , 
1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 

THREAT 
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HINCY=0.0 

RETURN 

END 

C 

FUNCTION HINCZ( I , J ,K) 

C 

C0M110N/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 

COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,XO(29) ,YO(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZOI(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY I , NZ 1 , N , M , MQ , DT , XMU , E PSO , EPS , NPTLIM , 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX.TY.TZ, AMP, ALPHA, BETA, IDLS 

THREAT 

HINCZ=0 .0 
RETURN 
END 

FUNCTION EINCXP(XP,YP,ZP) 

COMMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 

COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,ZO(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX, NY, NZ ,NX1 , NY1 ,NZI ,N,M,MQ,DT, XMU, EPSO, EPS, NPTLIM, 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 

THREAT 

EINCXP =0.0 
IF(IDLS.EQ.l) GO TO 1 
TAU = T - (UY0(23) - YP)/C 
IF (TAU.LE.O) GO TO 1 

EINCXP = AMP * ( EXP ( -ALPHA* TAU ) - EXP(-BETA*TAU) ) 

1 CONTINUE 
RETURN 
END 
C 

FUNCTION EINCYP(XP,YP,ZP) 

C 

COMMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 

C0M1T0N/ GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,YO(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX, NY, NZ ,NX1 ,NY1 ,NZI ,N,M,MQ,DT, XMU, EPSO, EPS, NPTLIM, 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPFAC, IP, TX.TY.TZ, AMP, ALPHA, BETA, IDLS 

C 

EINCYP = 0.0 
RETURN 
END 
C 

FUNCTION EINCZP(XP, YP,ZP) 

C 

COMMOH/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZO(29) 
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COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,ZO(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DXO(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0l(28) ,DY0I(28) ,DZOI(28) 

COMMON/ EXTRAS/NX » NY, NZ ,NXl ,NY1 ,NZl , N,M,MQ,DT, XMU, EPSO, EPS, NPTLIM, 
1 NN.NPTS.LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
E1NCZP = 0.0 
RETURN 
END 
C 

FUNCTION HINCXP(XP, YP,ZP) 

C 

COMMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 
C0MM0N/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,YO(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX, NY, NZ ,NX1 ,NY1 ,NZI ,N,M,MQ,DT,XMU,EPSO, EPS, NPTLIM, 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPF AC , IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
C 

HINCXP =0.0 
RETURN 
END 
C 

FUNCTION HINCYP(XP, YP , ZP ) 

C 

COMMON/ UGRID/UX( 28) ,UY(28) ,UZ(28) ,UX0(29) ,UY0(29) ,UZ0(29) 
C0MM0N/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI (29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY1 , NZ 1 , N , M, MQ , DT , XMU , EPSO , EPS , NPTLIM , 
1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPF AC , IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
C 

1IINCYP = 0.0 
RETURN 
END 
C 

FUNCTION UINCZP(XP , YP, ZP) 

C 

COMMON/UGRID/UX(28) ,UY(28) ,UZ(28) ,UX0(29) ,UYO(29) ,UZ0(29) 
C0MM0N/GRID/X(28) ,Y(28) ,Z(28) ,X0(29) ,Y0(29) ,Z0(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DX0(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZ0I(28) 
COMMON/EXTRAS/NX, NY, NZ ,NX1 ,NY1 ,NZ1 ,N,M,MQ,DT, XMU, EPSO, EPS, NPTLIM, 

1 NN, NPTS, LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX.TY.TZ, AMP, ALPHA, BETA, IDLS 
C 

THREAT 

UINCZP=0 »0 
RETURN 
END 
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3.5.8 Subroutine SAVESB 


The subroutine SAVESB saves the tangential E-field components on the 
subboundary selected for the M=1 loop, when the model is unexpanded. The 
subboundary is defined by the values set in a data statement for INEAR, JNEAR , 
KNEAR, IFAR, JFAR, KFAR. The values must be consistent with the expansion 
factor value selected. For example, an expansion factor of 4 (2,4,7 and 14 are 
allowed for a problem space 28 cells on a side) require seven cells on a side 
in the subboundary volume. In this case the values of INEAR and IFAR would be 
N and N+7, where N is an integer from 1 to 22. 

SAVESB treats opposing sides of the subboundary in succession (INEAR 
and IFAR, JNEAR and JFAR, KNEAR and KFAR sides respectively), saving for each 

side the two E-field components tangential to the surface (EY and EZ for the 

INEAR and IFAR sides for example) in the ARAY (L, IT). These components on a 

seven cell on a side subboundary form 9x8 and 8x9 arrays on a side. This 

provides field components on the edges of the subboundary (8 field components 
field components on each cell edge) or bracketing the edges of the subboundary 
(9 field components across and one cell on either side of the face at each cell 
center). As a result, in OUTBND it is possible to interpolate field values 
between the ARAY stored values using a single simple interpolation scheme. If 
the components running 9 across and about a face were not saved, rather only 7 
across the face, then an extrapolation scheme for the fields along the edges 
would be required, unnecessarily complicating matters. 

With 9x8 and 8x9 arrays of two field components on a face and 6 faces, 
the total number of stored components for a single call to SAVESB is 864, when 
an expansion factor of 4 is selected. SAVESB is called for each time step, 
which is labeled by IT in SAVESB, and ARAY is dimensioned as ARAY (864,500), 
when the expansion factor is 4 and no more than 500 time steps are required in 
the M=1 loop. These dimensions must be manually reset for other cases. 

Figure 13 is a flow chart of SAVESB. 
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* L Is advanced by one for every stored value 

Figure 12. SAVESB Flow Chart 
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SUBROUTINE SAVESB 
C 

COMMON/EFIELD/EX(29,29,29) ,EY(29,29,29) ,EZ(29, 29,29) 

COMMON/ EXTRAS/ NX , NY , NZ , NX1 , NY 1 , NZ 1 , N , M , MQ , DT , XMU , EPS 0 , EPS , NPTLIM , 
1 NN , NPTS , LMAX , SIGMA, C , T , PI , EXPF AC , I P , TX , TY , TZ , AMP , ALPHA , BETA , IDLS 
C 

COMMON/ RAY/ ARAY( 864 , 500) 

C 

DATA INEAR , JNEAR , KNEAR , IFAR, JFAR, KFAR/ 7,9,12,14,16,19/ 

THIS SUBROUTINE SAVES THE TANGENTIAL E FIELD COMPONENTS ON THE 
SUBBOUNDARY 


EXPFAC=4.0 
LIMIT=28/ EXPF AC 
L=1 

IF(N.NE.O) GO TO 10 

IT=1 

PRINT 5 

5 FORMAT(*SAVESB CALLED*) 
10 CONTINUE 

IMIN=INEAR 

JMIN=JNEAR 

K1-1IN=KNEAR 

EY ON INEAR, J,K SIDE 

IMAX-IMINfLIMIT 
JMAX= Jt II N+L IMIT 
KMAX=KMIN+LIMIT 

IL=IMIN— 1 
JL=JMIN-1 
KL-KMIN-1 

DO 90 J=JL, JMAX 
DO 90 K=KMIN,KMAX 
ARAY(L,IT)=EY(INEAR,J,K) 
L=L+1 

90 CONTINUE 

EY ON IFAR, J,K SIDE 

DO 110 J=JL, JMAX 
DO 110 K=KMIN,KMAX 
ARAY(L,IT)=EY(IFAR,J,K) 
L=L+1 

110 CONTINUE 

EZ ON INEAR, J,K SIDE 


DO 130 J=JMIN, JMAX 
DO 130 K=KL , KMAX 
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non oooo oon noon ooo 


AKAY(L,IT)=EZ( INEAR, J,K) 
L=L+1 

130 CONTINUE 

EZ ON IFAR, J ,K SIDE 

DO 150 J-JMIN, JMAX 
DO 150 K-KL.KMAX 
ARAY(L,IT)=EZ( IFAR, J ,K) 
L=L+1 

150 CONTINUE 

EX ON I , J , KNEAR SIDE 


DO 170 I=IL, IMAX 
DO 170 J=JMIN, JMAX 
ARAY(L,IT)-EX(I,J, KNEAR) 
L=L+1 

170 CONTINUE 

EX ON I , J , KF AR SIDE 

DO 190 I=IL , IMAX 
DO 190 J=JMIN, JMAX 
ARAY(L, IT)=EX( I, J ,KFAR) 
L=L+1 

190 CONTINUE 

EY ON I, J, KNEAR SIDE 


DO 210 I=IMIN, IMAX 
DO 210 J=JL, JMAX 
ARAY ( L , IT)=EY( I , J , KNEAR) 
L=L+1 

210 CONTINUE 

EY ON I , J , KF AR SIDE 

DO 230 I=IMIN, IMAX 
DO 230 J=JL, JMAX 
ARAY(L, IT)=EY( I, J ,KFAR) 
L=L+1 

230 CONTINUE 
C 

C EX ON I, JNEAR.K SIDE 

C 

C 

DO 250 I=IL, IMAX 
DO 250 K=KMIN,KMAX 
ARAY ( L , IT ) =EX( I , JNE AR , K) 
L=L+1 

250 CONTINUE 
C 
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non noon on 


EX ON I,JFAR,K SIDE 

DO 270 I=IL , IMAX 
DO 270 K=KMIN,KMAX 
ARAY(L, IT)=EX( I , JFAR,K) 
L=L+1 

270 CONTINUE 

EZ ON I,JNEAR,K SIDE 


DO 290 I=IMIN, IMAX 
DO 290 K=KL, KMAX 
ARAY(L,IT)=EZ(I, JNEAR,K) 
L=L+1 

290 CONTINUE 

C EZ ON I, JFAR,K SIDE 

DO 310 I=IMIN, IMAX 
DO 310 K-KL.K11AX 
ARAY(L,IT)=EZ(I, JFAR,K) 
L=L+1 

IF(IT.NE.l) GO TO 310 
310 CONTINUE 
C 

IT-IT+1 

RETURN 

END 
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3.5.9 


Subroutine OUTBND 


Subroutine OUTBND in the M=2 loop reads the data from ARAY acquired in 
SAVESB in the M=1 loop via a common statement. The data is first interpolated 
in time with each time step replaced by an expansion factor (EXPFAC or as in 
the present code i|) of time steps. Next the data is interpolated spatially 
between the unexpanded grid points (UX, UY, UZ, UXO, UYO, UZO of 
COMMON/UGRID/ . . . ) with a simple linear scheme and the expanded grid point 
locations (X, Y,Z ,XO,YO,ZO of COMMON /GR I D/. ..) . The interpolated data is placed 
on the faces of the problem space in the M=2 loop and acts in lieu of the 
radiation boundary condition defined in ERAD that is used on the M=1 loop. The 
same sequence of faces and components as used in SAVESB is used in OUTBND to 
insure that the interpolated data appears on the correct face in the correct 
location. Figure Hi is a flow chart of OUTBND. 
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Figure 13. OUTBND Flow Chart 
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SUBROUTINE OUTBND 


C 

C THIS SUBROUTINE ESTABLISHES THE TANGENTIAL E-FIELD COMPONENTS 

C ALONG THE OUTER BOUNDARY OF THE SUBBOUNDARY 

C 

COMMON/EFIELD/EX(29,29,29),EY(29,29,29),EZ(29,29,29) 
COMMON/GRID/X(28) ,Y(28) ,Z(28) ,XO(29) ,YO(29) ,ZO(29) , 

1 DX(29) ,DY(29),DZ(29) ,DXO(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DYOI(28) ,DZ0I(28) 
COMMON/UGRID/UX( 28) ,UY(28) ,UZ( 28) ,UXO(29) ,UYO(29) ,UZ0(29) 

COMMON/ EXTRAS/ NX, NY, NZ.NX1 ,NY1 ,NZ1 ,N,M,MQ,DT,XMU,EPSO, EPS, NPTLIM, 

1 NN, NPTS,LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
C 

COMMON/ RAY / ARAY ( 864 ,500) 

DIMENSION ARRAY (864) , CRAY (9, 9) 

C 

DATA INEAR , JNEAR, KNEAR, IFAR, JFAR, KFAR/7 , 9,12,14,16,19/ 

IF(N.NE.O) GO TO 10 
PRINT 5 

5 FORMAT(*OUTBND CALLED*) 

IT=1 

10 CONTINUE 
E=M0D(N,4) 

C 

DO 50 LL=1 , 864 

ARRAY(LL)=ARAY(LL, IT)+(E/EXPFAC)*(ARAY(LL,IT+1)-ARAY(LL,IT) ) 

50 CONTINUE 
C 

IF(E.NE.3) GO TO 20 
IT=IT+1 
20 CONTINUE 
C 

LIMIT-28/ EXPF AC 
LONE-LIMIT+1 
LTW0=LIMIT+2 
LL=1 

DO 80 J=1 ,LTWO 
DO 80 K-l.LONE 
CRAY (J,K) -ARRAY (LL) 

LL-LL+1 
80 CONTINUE 
C 

DO 90 J-1,28 
JJ=INT((J+1)/EXPFAC)+1 
DO 90 K-1,28 
KK=INT((K-2)/EXPFAC)+l 
JS =J J+J NEAR- 1 
KS-KK+KNEAR- 1 

EY1=CRAY( JJ+1 ,KK)+( ( ZO(K)-UZO(KS) )/ (UZ0(KS+1)-UZ0(KS) ) ) 

1 *(CRAY(JJ+1 ,KK+1) -CRAY (JJ+1 ,KK)) 

EY2-CRAY( JJ, KK)+( ( Z0( K)-UZO( KS ) ) / ( UZ0( KS+1 )-UZO( KS ) ) ) 

1 *(CRAY ( J J ,KK+1 )-CRAY( JJ ,KK) ) 

EY( 1 , J,K)=EY2+( (Y( J)-UY( JS-1) )/ (UY( JS )-UY( JS-1) )) 

1 *(EY1-EY2) 

90 CONTINUE 
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DO 100 J=l,LTWO 
DO 100 K-l.LONE 
CRAY ( J , K)=ARRAY (LL ) 

LL=LL+1 
100 CONTINUE 
C 

DO 110 J=1 , 28 
JJ=INT((J+1)/EXPFAC)+1 

DO 110 K=1 ,28 
KK=INT((K-2)/EXPFAC)+l 
J S = J J+J NEAR- 1 
KS=KK+KNEAR-1 

EY1=CRAY( JJ+1 ,KK)+( (ZO(K)-UZO(KS) )/ (UZ0(KS+1)-UZ0(KS) ) ) 
1 *( CRAY( JJ+1 ,KK+1 )-CRAY( JJ+1 ,KK) ) 

EY2=CRAY( JJ ,KK)+( (ZO(K)-UZO(KS) )/(UZO(KS+l)~UZO(KS) ) ) 

1 * ( CRAY ( J J , KK+1 ) -CRAY ( J J , KK) ) 

EY(29 , J,K)=EY2+( (Y( J )-UY( JS-1) )/(UY( JS)-UY( JS-1) ) ) 

1 *(EY1-EY2) 

110 CONTINUE 
C 

DO 120 J=1 , LONE 
DO 120 K=1 ,LTUO 
CRAY( J ,K)=ARRAY(LL) 

LL=LL+1 
120 CONTINUE 
C 

DO 130 J=1 , 28 

JJ=INT( ( J-2)/EXPFAC)+l 

DO 130 K-1,28 

KK=INT((K+1)/EXPFAC)+1 

JS=JJ+JNEAR-1 

KS=KK+KNEAR-1 

EZ1=CRAY( JJ+1 ,KK)+( (Z(K)-UZ(KS-l) )/ (UZ(KS)-UZ(KS-l) ) ) 

1 * (CRAY( JJ+1, KK+1 )-CRAY( JJ+1, KK)) 

EZ2=CRAY(JJ,KK)+( (Z(K)-UZ(KS-l ) )/ (UZ(KS)-UZ(KS-l) )) 

1 *(CRAY( JJ ,KK+1 )-CRAY( JJ ,KK) ) 

EZ( 1 , J ,K)*»EZ2+( (Y0( J)-UYO( JS) )/ (UYO( JS+1)-UY0( JS) ) ) 

1 *(EZ1-EZ2) 

130 CONTINUE 
C 

DO 140 J=1 , LONE 
DO 140 K=1 ,LTWO 
CRAY ( J , K)=ARRAY( LL) 

LL-LL+1 
140 CONTINUE 
C 

DO 150 J=1 , 28 
JJ=INT((J-2)/EXPFAC)+l 
DO 150 K= 1 , 2 8 
KK=INT((K+1)/EXPFAC)+1 
JS=J J+JNEAR-1 
KS =KK+KNEAR- 1 

EZ1=CRAY (JJ+1 ,KK)+( (Z(K)-UZ(KS-l ) )/ (UZ(KS)-UZ(KS— 1 ) ) ) 

1 * ( CRAY (JJ+1, KK+1 ) -CRAY ( J J+l , KK) ) 
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EZ2=CRAY( JJ ,KK)+( (Z(K)-UZ(KS-l) )/ (UZ(KS)-UZ(KS-l ) ) ) 

1 * ( CRAY ( J J , KK+1 ) -CRAY ( J J , KK) ) 

EZ(29,J,K)=EZ2+((Y0(J)-UY0(JS))/(UY0(JS+1)-UY0(JS))) 

1 *(EZ1-EZ2) 

150 CONTINUE 
C 

DO 160 1=1 , LTWO 
DO 160 J=1,L0NE 
CRAY(I,J)=ARRAY(LL) 

LL=LL+1 
160 CONTINUE 
C 

DO 170 1=1,28 
II=INT( (1+1) /EXPFAC)+1 
DO 170 J=1 ,28 
JJ=INT( ( J-2) /EXPFAC)+1 
IS=I I+INEAR-1 
J S = J J+J NE AR- 1 

EX1=CRAY(II , JJ+1 )+( (X(I)-UX( IS-1) )/ (UX(IS)-UX( IS-1) ) ) 
1 *(CRAY(II+1 , JJ+1 )-CRAY ( II , JJ+1 ) ) 

EX2=CRAY(II,JJ)+((X(I)-UX(IS-1))/(UX(IS)-UX(IS-1))) 

1 *(CRAY(II+1,JJ)-CRAY(II,JJ)) 

EX(I , J , l)=EX2+( (Y0( J)-UYO( JS) )/ (UYO( JS+1 )-UYO( JS) ) ) 

1 *(EX1-EX2) 

170 CONTINUE 
C 

DO 180 1=1, LTWO 
DO 180 J=1,L0NE 
CRAY ( I , J ) = ARRAY ( LL ) 

LL=LL+1 
180 CONTINUE 
C 

DO 190 1=1,28 
II=INT( ( 1+1) / EXPFAC)+1 
DO 190 J=1 ,28 
JJ=INT( ( J-2 ) /EXPFAC)+1 
IS=I I+INEAR-1 
J S = J J+ J NEAR- 1 

EX1=CRAY( II , JJ+1 )+( (X( I)-UX( IS-1 ) )/ (UX(IS)-UX(IS-l) )) 
1 * ( CRAY ( I 1+1 , J J+l ) -CRAY (II, JJ+1)) 

EX2=CRAY( II, JJ)+((X(I)-UX( IS-1 ))/(UX(IS)-UX( IS-1))) 

1 *(CRAY(II+1 , JJ)-CRAY(II,JJ)) 

EX(I , J , 29)=EX2+( (Y0( J)-UYO( JS) )/ (UYO( JS+1 )-UYO( JS) ) ) 

1 *(EX1-EX2) 

190 CONTINUE 
C 

DO 200 1=1, LONE 
DO 200 J=1,LTW0 
CRAY(I , J)=ARRAY(LL) 

LL=LL+1 
200 CONTINUE 
C 

DO 210 1=1,28 
II=INT( ( 1-2) /EXPFAC)+1 
DO 210 J=1 ,28 


110 



JJ=INT( ( J+l ) /EXPFAO+l 
IS=I I+INEAR-1 
JS=JJ+JNEAR-1 

EY1=CRAY(II, JJ+1 )+( (X0( I)-UX0( IS) )/ (UX0(IS+1)-UX0(IS) ) ) 
1 *( CRAY( II+l , JJ+1 )-CRAY(II, JJ+1)) 

EY2=CRAY(II, JJ)+( (XO(I)-UXO(IS) )/ (UX0(IS+1)-UX0(IS) )) 

1 *(CRAY(II+1 , JJ)-CRAY( II, JJ) ) 

EY( I , J , l)=EY2+( (Y( J)-UY( JS— 1 ) )/ (UY( JS)-UY( JS— 1) ) ) 

1 *(EY1-EY2) 

210 CONTINUE 
C 

DO 220 1=1, LONE 
DO 220 J=1 ,LTWO 
CRAY( I , J)=ARRAY(LL) 

LL=LL+1 
220 CONTINUE 
C 

DO 230 1=1,28 
II=INT((I-2)/EXPFAC)+l 
DO 230 J=1 , 28 
JJ=INT( ( J+l )/EXPFAC)+l 
IS=I I+INEAR-1 
JS=JJ+JNEAR-1 

EY1=CRAY( II , JJ+1 )+( (XO(I)-UXO( IS) )/ (UXO( IS+1 )-UXO( IS) ) ) 
1 *( CRAY( II+l, J J+l )-CRAY( II, JJ+1)) 

EY2=CRAY( II , JJ)+( (XO(I)-UXO( IS) )/ (UX0(IS+1)-UX0( IS) ) ) 

1 *( CRAY( II+l , JJ)-CRAY(II, JJ) ) 

EY( I , J,29)=EY2+( (Y(J )-UY( JS-1) )/(UY( JS)-UY( JS-1) )) 

1 *(EY1-EY2) 

230 CONTINUE 
C 

DO 240 1=1 ,LTWO 
DO 240 K=1,L0NE 
CRAY ( I , K )= ARRAY ( LL ) 

LL=LL+1 
240 CONTINUE 
C 

DO 250 1=1,28 
II=INT((I+1)/EXPFAC)+1 
DO 250 K= 1 , 2 8 
KK=INT((K-2)/EXPFAC)+l 
IS=I I+INEAR-1 
KS=KK+KNEAR- 1 

EX1=CRAY( II, KK+l)+( (X(I)-UX( IS-1 ) )/ (UX( IS)-UX(IS-l) ) ) 

1 * ( CRAY ( I 1+1 , KK+1 ) -CRAY ( I I , KK+1 ) ) 

EX2=CRAY( II , KK)+( (X( I)-UX( IS-1 ) )/ (UX( IS)-UX( IS-1 ) ) ) 

1 *(CRAY( II+l ,KK)-CRAY(II,KK) ) 

EX( I, 1 ,K)=EX2+( (ZO(K)-UZO(KS) )/ (UZ0(KS+1 )-UZO(KS) ) ) 

1 *(EX1-EX2) 

C 

250 CONTINUE 

DO 260 1=1 ,LTWO 
DO 260 K=1 ,LONE 
CRAY ( I , K )= ARRAY ( LL ) 

LL=LL+1 
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260 CONTINUE 


DO 270 1=1,28 
II=INT((I+1)/EXPFAC)+1 
DO 270 K=i ,28 
KK=INT((K-2)/EXPFAC)+l 
IS=II+INEAR-1 
KS=KK+KNEAR- 1 

EXl=CRAY(II,KK+l)+( (X( I)-UX(IS-1 ) )/ (UX(IS)-UX(IS-I) ) ) 

1* ( CRAY ( I 1+1 , KK+I ) -CRAY ( I I , KK+1 ) ) 

EX2=CRAY(II,KK)+((X(I)-UX(IS-1))/(UX(IS)-UX(IS-I))) 

I *(CRAY(II+1,KK)-CRAY(II,KK)) 

EX( I , 29 ,K)=EX2+( (ZO(K)-UZO(KS) )/ (UZ0(KS+1 )-UZO(KS) )) 

1 *(EX1-EX2) 

270 CONTINUE 
C 

DO 280 1=1, LONE 
DO 280 K-I.LTWO 
CRAY(I,K)=ARRAY(LL) 

LL=LL+1 
280 CONTINUE 
C 

DO 290 1=1,28 

II=INT( ( 1-2 ) /EXPFAC)+1 

DO 290 K=l,28 

KK=INT((K+1)/EXPFAC)+1 

IS=II+INEAR-1 

KS=KK+KNEAR-1 

EZl=CRAY(II,KK+l)+( (XO(I)-UXO( IS) )/ (UX0(IS+1)-UX0(IS) ) ) 
1*( CRAY(II+1 , KK+1 ) -CRAY ( I I, KK+1 ) ) 

EZ2=CRAY(II,KK)+( (XO(I)-UXO(IS) )/(UXO(IS+l )-UX0(IS) ) ) 

1 * ( CRAY ( I 1+1 , KK) -CRAY ( I I , KK) ) 

EZ( 1 , 1 , K)=EZ2+( (Z(K)-UZ(KS-l ) ) / ( UZ( KS)-UZ( KS-1 ) ) ) 

1 *(EZ1-EZ2) 

290 CONTINUE 
C 

DO 300 1=1, LONE 
DO 300 K=1 ,LTWO 
CRAY ( I , K ) =ARRAY ( LL ) 

LL=LL+1 
300 CONTINUE 
C 

DO 310 1=1,28 
II=INT((I-2)/EXPFAC)+l 
DO 310 K= 1 , 2 8 
KK=INT( (K+l ) /EXPFAC)+1 
IS=II+INEAR-1 
KS=KK+KNEAR- 1 

EZl=CRAY(II , KK+1 )+( (XO(I)-UXO(IS) )/ (UX0(IS+1)-UX0(IS) ) ) 
1*(CRAY(II+1, KK+1 )-CRAY( II, KK+1)) 

EZ2=CRAY(II ,KK)+( (X0( I)-UXO(IS) )/ (UX0(IS+1)-UX0(IS) ) ) 

1 *(CRAY(II+1 ,KK)~CRAY(II,KK)) 

EZ(I, 29 ,K)=EZ2+( (Z(K)-UZ(KS-l ) )/(UZ(KS)-UZ(KS-l) ) ) 

1 *(EZ1-EZ2) 

310 CONTINUE 
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RETURN 

END 


3.5.10 Subroutine DATASAV 


Subroutine DATASAV is called by the main program DRIVER for each loop 
(M=1 ,2) selected (by MM setting) as frequently as output values of the induced 
surface charges (n*E) and currents (nxH) at the selected test points are 
required (determined by IP). A flow chart for the M=1 loop is given to 
Figure 15. 

On the first call in a particular loop to Subroutine DATASAV 
parameters that describe the test points are set. These parameters are defined 
in Table 6. For this and all subsequent calls within the loop, DATASA performs 
the operations described in the following paragraphs for each of the three 
field components from which the charge and the two current components are 
calculated. These operations allow one to define an arbitrary test point 
location and orientation where current and charge responses can be found from 
interpolated and extrapolated field responses at the appropriate surrounding 
I,J,K locations of the calculated fields. 


Table 6. Variables Input in Subroutine DATASAV 


Variable 

NAME 

NPTS 

LOC(NPT) 


Description 

Number of test points 

fch 

Test point number for the M test point 


NPLANE (NPT ) 


XOBS(NPT) 
YOBS (NPT) 
ZOBS (NPT ) 
THETA (NPT) 


Scattering Plane indicates |Jie orientation 
of the plane on which the M c test point is 
located 

1 - Normal is in the +x direction 

2 - Normal is in the -x direction 

3 - Normal is in the +y direction 

4 - Normal is in the -y direction 

5 - Normal is in the +z direction 

6 - Normal is in the -z direction 

th 

x coordinate of the M test point 

til 

y coordinate of the M test point 

tfl 

z coordinate of the M test point 

Measurement orientation in degrees from 
one of the principle axes of the problem 
space 
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Figure 14. 


Flow Chart for Subroutine DATASAV 
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Using the test point location and the scattering plane orientation, 
eight field component locations are selected. These field component locations 
are distributed about the test point as illustrated below in Figure 1 6. 



Figure 15 . Distribution of Field Component Locations About Test Point 

The selection process compared all the field component locations 
against the test point location in one of the three dimensions using a 
criterion of the form 

if XOBS(NPT) < X(I+1), then I max = (1+1) and I fflln = (I) 

(here the x coordinate has been used for illustrative purposes only), and in 
the remaining dimensions using a criterion of the form 

if YOBS(NPT) < Y ( J+ 1 ) , the J , = (J) and J = (J+1). 

min max 

By letting the coordinates by X(I) or X q (I), Y(J) or Y q (J), and Z(K) or Z q (K) 
in the above criteria, depending on the field component selected, the test 
point is surrounded by the eight field components. 

The eight field component values associated with the eight field 
component locations are given by: 
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where 5 is any one of the six field components. From these field components 
the value of 5(X0BS(NPT), YOBS(NPT), ZOBS(NPT)) can be found by a combination 
of interpolation and extrapolation. The exact order in which the components 
are combined depends on the test point location and NPLANE. An example using 
the field component arrangement shown in Figure 17 is given to illustrate the 
process. 

First an interpolation is performed in the z direction to yield 


£(I . ,J . ,K . ) + (5(1 , ,J . ,K ) 

s min min mm min min max 

Az 

I 2 t ^ 0 t K . ) ) * * rj 

min min min AZ 

5(1 ,J . ,K . ) + (5(1 ,J . ,K ) 

s max mm mm max min’ max 

Az 

^ "'‘max’ J min’ ^min^ AZ 
^ ■''min’ ^max’ K min^ + ^ '*'min’ lI max' K max^ 


' « I . ln . J max» K « i n )) * AZ~ 

^4 _ ^ I max’ ^max' K min^ + ^ *max’^max’*Snax^ 

At 

- 5(1 ,J ,K . )) • ~ 
s max max min AZ 
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These interpolated values are distributed about the test point as shown in 
Figure 18. They are then interpolated in the x direction to yield 


12 

" *1 

+ ( ^2 

- V 

Ax 

AX 

34 


+ U 3 


Ax 

AX 


Finally, the values and as shown in Figure 19, are extrapolated to the 

test point location to give £(X0BS(NPT), YOBS(NPT), ZOBS(NPT)), using the 
expression 

£( XOBS( NPT) , YOBS(NPT), ZPBS(NPT)) = ^ . 



Figure 16. Example Field Component Arrangement about a Test Point 
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Figure 17. Interpolated Field Values about a Test Point 


To each value of the scattered field component at the test location, 
£(X0BS(NPT), YOBS(NPT), ZOBS(NPT)), there is added the corresponding incident 
field component at the test point location. This is accomplished by adding the 
appropriate function value of EINCXP, EINCYP, EINCZP, HINCXP or HINCZP. The 
resulting total field components correspond to the induced surface charges and 
currents. 



( 6 Test Point (XOBS(NPT) ,Y0BS(NPT) .ZOBS(NPT) 


Figure 18. Interpolated Field Values Used to 
Extrapolate to Test Point Location 
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Two current components are calculated on any surface of the scatterer. 
These components can be combined to yield other orthogonal current component 
pairs corresponding to those measuring in an experiment where the current 
sensors are not aligned in the direction of the x, y or z axis. This situation 
is commonly the case for measurements made on aircraft wings where the current 
sensors are aligned with the axis of the wing, which generally is canted with 
respect to one of the principal axes of the problem space. 

The angle THETA(NPT) is the sensor rotation at test point with respect 
to one of the principal axes. It is shown graphically in Figure 20. Simply 
stated, 0 is always defined as falling in the first quadrant of the plane in 
which the scattering surface falls. As a result of this definition, the two 
rotated, orthogonal, current components, called HI and H2 , are always given by: 

HI = Acos0 + Bsin© 

H2 = -Asin0 + BcosG 

where 

A = Hz, B = Hy for NPLANE = 1 or 2 

A = Hz, B = Hx for NPLANE = 3 or 4 

A = Hx, B = Hy for NPLANE = 5 or 6 

Those variables used locally in Subroutine DATASAV are defined in 

Table 7. 

front side 



Figure 19. Graphical Description of THETA 

(Angle defining sensor rotation). 
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Table 7. Local Variables - Subroutine DATASAV 


Variable 

Name Description 

A A principle H field component 


B 


A principle H field component 


COSAN 

EX1 ,EX2 
EX3.EX4 
EX1 2,EX34 

EY1 , EY2 
EY3.EY4 
EY1 2.EY34 

EZ1 , EZ2 
EZ3,EZ4 
EZ1 2, EZ34 

HX1 , HX2 
HX3.HX4 
HX1 2, HX3 4 

HY1 ,HY2, 

HY3,HY4 

HY12.HY34 

HZ1 ,HZ2 
HZ3.HZ4 
HZ1 2 , HZ3 4 


} 

} 

) 

} 

} 

} 


Cosine of THETA 

Interpolated E^ field components 

Interpolated E^ field components 

Interpolated E field components 
z 

Interpolated field components 

Interpolated field components 

Interpolated H field components 
z 


L 

LOC 

NPLANE 

NPTS 

IM, IX n 
JM, JX > 
KM.KX J 


I MIN , IMAX n 
JMIN, JMAX \ 
KMIN.KMAX J 

SINAN 


TERM1 , TERM2 
TERM3 , TERM4 
TERM5 , TERM6 
TERM7, TERM8 
TERM9 


Counter on the number of time steps output 

Test point number (array dimensioned on the 
number of test points) 

Scattering Plane (array dimensioned on the 
number of test points) 

Number of test points 


Values of I . ,1 ,J . ,J ,K . ,K 

min max min max nun max 


Same as above (arrays dimensioned on the 
number of test points) 

Sine of THETA 


Interpolation and extrapolation terms 
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Table 7. (continued) 


i 


Variable 

Name 

Description 

THETA 

Angle (degrees) of current sensor rotation 
(array dimensional on the number of test points) 

XBS 'N 

YBS > 
ZBS ' 

Coordinates of test point location 

XOBS > 

YOBS > 

Same as above (arrays dimensioned on the 

ZOBS ' 

number of test points) 
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SUBROUTINE DATASAV 
C 

COMMON/ EFIELD / EX(29,29,29) , EY(29,29,29) , EZ(29, 29,29) 

COMMON/ HFIELD / liX(29 ,29 ,29) , HY(29,29,29) , IiZ(29,29,29) 

COMMON/ EXTRAS/ NX, NY , NZ ,NX1 , NYI ,NZ 1 ,N , M ,MQ, DT, XMU, EPSO,EPS , NPTLIM, 
1 NN ,NPTS ,LMAX, SIGMA, C ,T , P I .EXPFAC , IP ,TX,TY,TZ , AMP, ALPHA, BETA, IDLS 
COMMON/ GRID/X(28) ,Y(28) ,Z(28) ,XO(29) ,YO(29) ,ZO(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DXO(28) ,DYO(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZ1(29) ,DX0I(28) ,DYOI(28) ,DZ0I(28) 

COMMON/ OUT/ESTORE(250 ,24) ,HSTOR1(2500,24) ,USTOR2(250,24) , 

1 TST0RE(250) 

COMMON/ CUMMUL/CUMM( 24) ,LOC(24) 

C 

DIMENSION X0BS(30) ,Y0BS(30) ,Z0BS(30) 

DIMENSION SDX0(28) ,SDYO(28) ,SDZ0(28) 

DIMENSION IMIN(30) , IMAX(30) ,JMIN(30) ,JMAX(30) ,KMIN(30) ,KMAX(30) 
DIMENSION THETA(30) ,NPLANE(30) 

DIMENSION SXO( 28) ,SY0(28) ,SZO(28) ,SX(28) ,SY(28) ,SZ(28) 

DIMENSION SEX(28,28,28) ,SEY(28,28,28) ,SEZ(28,28,28) 

C 

PIF=3. 1415926336/ 180.0 
L=N/IP 
TSTORE(L)=T 
C 

IF(L.NE.l) GO TO 1120 

IF(M.NE.l) GO TO 200 

LOC(l)=l 

LOC(2)=2 

LOC(3)=3 

L0C(4)=4 

LOC(5)=5 

LOC(6)=6 

NPLANE( 1)=3 

NPLANE(2)=4 

NPLANE(3)=1 

NPLANE(4)=1 

NPLANE(5)=I 

NPLANE ( 6) = 1 

XOBS( I)=-I . 8 

YOBS( 1)=0. 0 

ZOBS( 1 )=0 . 0 

XOBS(2)=-4.8 

YOBS(2)=-l • 92 

ZOBS(2)=0.0 

XOBS(3)=-l . 95 

YOBS(3)=-l . 68 

Z0BS(3)=0.42 

XOBS(4)=-l . 95 

YOBS(4)=-1.62 

Z0BS(4)=0.48 

X0BS(5)=-1 .95 

YOBS(5)=-I.68 

ZOBS(5)=0.54 

XOBS(6)=-I . 95 

YOBS(6)=-1.74 


123 



non 


ZOBS(6)=0.48 
THETA(1)=0.0 
THETA(2)=0.0 
TilETA(3)=90.0 
THETA(4)=0.0 
THETA(5)=90.0 
THETA(6)=0.0 
NPT3=6 
GO TO 499 
C 

200 CONTINUE 

IF(M.NE.2)GO TO 300 
L0C( 1)=7 
LOC(2)=8 
LOC(3)=9 
L0C(4)=10 
L0C(5) = U 
LOC(6)=I2 
NPLANE(1)=3 
NPLANE(2)=4 
NPLANE(3)=1 
NPLANE(4)=1 
NPLANE(5)=I 
NPLANE(6)=I 
X0BS(I)=-1.8 
YOBS(1)=0.0 
Z0BS( 1 )=0. 0 
X0BS(2)=-4.8 
Y0BS(2)=-1 . 92 
Z0BS(2)=0. 0 
XOBS(3)=-I.95 
YOBS(3)=-l .68 
ZOBS(3)=0.42 
XGBS(4)=-1 . 95 
Y0BS(4)=-1.62 
ZOBS(4)=0.48 
XOBS(5)=-1.95 
YOBS(5)=-1.68 
Z0BS(5)=0.54 
XOBS(6)=-1.95 
Y0BS(6)=-1 . 7 4 
Z0BS(6)=0.48 
THETA(1)=0.0 
THETA(2)=0.0 
THETA(3)=90.0 
THETA(4)=0.0 
THETA(5)=90.0 
THETA(6)=0.0 
NPTS=6 
GO TO 499 
C 

300 CONTINUE 

DETERMINE CELLS FOR EXTRAPOLATION 
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c 


1000 

1020 


1025 


1040 

1060 


1070 

1080 

1100 

1110 


C 

C 

C 

C 

c 

c 

G 

C 


1120 


DO 1110 NPT-l.NPTS 

IF(NPLANE(NPT) .GE.3) GO TO 1025 
DO 1000 1=1, NX1 

IF(X0B5(NPT) . LT.X(I+1) ) GO TO 1020 

CONTINUE 

CONTINUE 

IF(NPLANE(NPT).EQ.l) IMIN(NPT) = 1+1 
IF(NPLANE(NPT) .EQ.2) IMIN(NPT) = 1-1 
IMAX(NPT) = IMIN(NPT) + 1 
GO TO 1110 
CONTINUE 

IF(NPLANE(NPT) .GT.4) GO TO 1070 
DO 1040 J=1 , NY1 

IF(YOBS(NPT) .LT.Y(J+1) ) GO TO 1060 

CONTINUE 

CONTINUE 

IF(NPLANE(NPT) .EQ.3) JHIN(NPT) = J+l 

IF(NPLANE(NPT ) .EQ.4) JllIN(NPT) = J-l 

JMAX(NPT) = JMIN(NPT) + 1 

GO TO 1110 

CONTINUE 

DO 1080 K=1 , NZ1 

IF(ZOBS(NPT) . LT.Z(K+1) ) GO TO 1100 
CONTINUE 

IF ( NPLANE ( N PT ) . EQ . 5 ) KMIN(NPT) = K+l 
IF(NPLANE(NPT).EQ.6) KMIN(NPT) = K-l 
KMAX(NPT) = KMIN(NPT) + 1 
CONTINUE 

PRINT 407 , (XOBS(NPT) ,NPT=1 ,NPTS) 
PRINT 408 , (YOBS(NPT) ,NPT=1,NPTS) 
PRINT 409, (ZOBS(NPT) ,NPT=1 ,NPTS) 


NPLANE 

=2 

1 

IS 

THE 

YZ 

PLANE , 

NORMAL 

NPLANE 

= 

2 

IS 

THE 

YZ 

PLANE , 

NORMAL 

NPLANE 

= 

3 

IS 

THE 

XZ 

PLANE, 

NORMAL 

NPLANE 

= 

4 

IS 

THE 

XZ 

PLANE, 

NORMAL 

NPLANE 

= 

5 

IS 

THE 

XY 

PLANE, 

NORMAL 

NPLANE 

K 

6 

IS 

THE 

XY 

PLANE , 

NORMAL 


CONTINUE 

DO 1400 NPT=1 ,NPTS 
IM = IMIN(NPT) 

IX = IMAX(NPT) 

JH = JMIN(NPT) 

JX = JMAX(NPT) 

KM = KMIN(NPT) 

KX = KMAX(NPT) 

XBS » XOBS(NPT) 


AFT (TOWARD INCREASING I) 

F ORMARD ( T OW ARD DECREASING I) 

UP (TOWARD INCREASING J) 

DOWN (TOWARD DECREASING J) 

RIGHT (TOWARD INCREASING K) 
NEGATIVE (TOWARD DECREASING K) 
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YBS = YOBS(NPT) 

ZBS = ZOBS(NPT) 

ANGLE = THETA( NPT ) *PIF 


IF(NPLANE(NPT) .GT.4) GO TO 1270 
IF(NPLANE(NPT) .GT.2) GO TO 1155 

NPLANE = 1 OR 2 -- YZ PLANES 

COMPUTE TERMS FOR CHARGE 


DETERMINE JM,JX,KM,KX FOR EX 
DO 1125 J=1 ,NY1 

IF(YOBS(NPT) .LT.Y0(J+1) ) GO TO 1128 
1125 CONTINUE 
1128 JM = J 

JX = J+l 
DO 1130 K=1 ,NZ1 

IF(ZOBS(NPT) .LT. Z0(K+1 ) ) GO TO 1131 

1130 CONTINUE 

1131 KM = K 
KX = K+l 

IF(L.EQ.IOO) PRINT 410, IM,IX, JM,JX,KM,KX 
TERM1 = (ZBS - ZO(KM) ) /DZ0( KM) 

TERI'11 IS SCALE FACTOR FOR EX(I,J,KX) 

TERM2 = (YBS - Y0( JM))/DY0( JM) 

TERM2 IS SCALE FACTOR FOR EX(I,JX,K) 

TERM3 = (XBS - X0( IM) ) /DX0( IM) 

TERM4 = 0.5*(DX0(IX)+DX0(IM) ) 

TERM5 = 0. 5* (DY0( JX)+DY0( JM) ) 

TERM6 = 0. 5*(DZ0(KX)+DZ0(KM) ) 

EX1=EX( 111 , JM , KM)+( EX( IM , JM , KX)-EX( IM , JM , KM) ) *TERM1 
EX2=EX( IM, JX, Kll)+( EX( IM, JX,KX)-EX( IM, JX, KM) ) *TERM1 
EX3=EX( IX , JM , Ki-l)+( EX( IX , JM , KX)-EX( IX , JM , KM) ) *TERM1 
EX4=EX( IX,JX,KM)+(EX( IX, JX,KX)-EX( IX, JX.Ki'l) )*TERM1 
EX12 = EX1 + (EX2-EX1)*TERM2 
EX34 = EX3 + ( EX4-EX3 ) *TERM2 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1134 

PRINT 514, ( ( (EX( I, J,K) , I=IM, IX) , J=JM, JX) , K=KM,KX) 
PRINT 515, EX1,EX2,EX3,EX4,EX12,EX34 
1134 CONTINUE 
C 

C COMPUTE TERMS FOR CURRENT DENSITY 

C 

C 

C DETERMINE KM,KX FOR HY, JM,JX SAME AS FOR EX 

C 

DO 1135 K=1 ,NZ1 
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IF(Z0BS(NPT).LT.Z(K+1)) GO TO 1136 

1135 CONTINUE 

1136 KM * K 
KX = K+l 

IF(L.EQ.IOO) PRINT 410, IM , IX, JM, JX, KM,KX 
TERM1 = (ZBS - Z0(KM))/DZ0( KM) 

TERM2 = (YBS - Y0( JM) ) /DY0( JM) 

C TERM2 IS SCALE FACTOR FOR I1Y(I,JX,K) 

TERM3 = (XBS - X0(IM))/DX0(IM) 

TERM4 = 0 . 5* ( DX0( IX)+DX0( IM) ) 

TERM5 = 0. 5*(DY0( JX)+DY0( JM) ) 

TERM6 = 0.5*(DZ0(KX)+DZ0(KM) ) 

TERM7 = (ZBS-Z(KIl) ) /TERM6 
C TERM7 IS SCALE FACTOR FOR HY(I,J,KX) 

HY1=HY( IM , JM , KM)+( UY( IM , JM , KX) -11Y( IM , JM , KM) ) *TEKM7 
liY2=HY(IM,JX,KM)+(UY(IM,JX,KX)-HY(IM,JX,KM))*TERM7 
11Y3=HY(IX, JM,KM)+(HY( IX, JM,KX)-1IY( IX, JM,KM) )*TERM7 
IlY4=UY(IX,JX,KM)+(HY(IX,JX,KX)-HY(IX,JX,Ki'l))*TERM7 
C 

HY12 = HY1 + (HY2-11Y1)*TERM2 
HY34 = HY3 + (HY4-HY3)*TERM2 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1137 

PRINT 511,(((HY(I,J,K),I=IM,IX),J=JM,JX),K-KM,KX) 
PRINT 512, HY1 , HY2 , HY3 , HY4 , HY1 2 , HY34 

1137 CONTINUE 


DETERMINE JM,JX,KM,KX FOR HZ 
DO 1138 J=1 ,NY1 

IF(YOBS(NPT) .LT.Y( J+l) ) GO TO 1140 
1138 CONTINUE 
1140 JM = J 
JX = J+l 
DO 1144 K= 1 , NZ 1 

IF( ZOBS(NPT) .LT. Z0(K+1 ) ) GO TO 1146 
1144 CONTINUE 
1146 KM = K 
KX = K+l 

IF(L.EQ.IOO) PRINT 410,IM,IX, JM, JX,KM,KX 
TERM1 = (ZBS - Z0(K11))/DZ0(KM) 

C TERM1 IS SCALE FACTOR FOR I1Z(I,J,KX) 

TERM2 = (YBS - YO( JM) )/DY0( JM) 

TERM3 = (XBS - X0(IM))/DX0(IM) 

TERM4 = 0. 5*(DX0( IX)+DXO( IM) ) 

TERM5 = 0. 5*(DY0( JX)+DYO(JM) ) 

TERM6 = 0.5*(DZ0(KX)+DZ0(KM) ) 

C 

TERM8 = (YBS-Y(JM) )/TERM5 
C TERMS IS SCALE FACTOR FOR HZ(I,JX,K) 

HZ 1=HZ ( IM , JM , KM)+( HZ( IM , JM , KX)-HZ( IM , JM , KM) ) *TERfll 
HZ2=HZ( IM, JX,KM)+(HZ( IM, JX,KX)-HZ( IM, JX,KM) )*TERM1 
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HZ3=HZ( IX, JM,KM)+( HZ( IX, JM,KX)-HZ( IX, JM,KM) )*TERM1 
HZ4=HZ(IX,JX,KM)+(HZ(1X,JX,KX)-HZ(IX,JX,KM))*TERM1 
HZ 12 = HZ1 + (HZ2-L1Z1 )*TERM8 
HZ34 = HZ3 + ( HZ4-HZ3 ) *TER118 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1148 

PRINT 504 , ( ( (HZ( I , J ,K) ,I=IM,IX) ,J=JM,JX) ,K»KM,KX) 
PRINT 505, HZ1,HZ2,HZ3,HZ4,HZ12,HZ34 
1148 CONTINUE 

IF(NPLANE(NPT) .EQ.2) GO TO 1150 

TERM9 = (X(IM)-XBS)/TERM4 

TER119 IS -SCALE FACTOR FOR (EX,HY,UZ) (IX, J ,K) 


EST0RE(L,NPT)=EX12 - (EX34-EX12)*TERII9 + EINCXP(XBS, YBS.ZBS) 
A=11Z12-(HZ34-HZ12)*TERM9+HINCZP(XBS , YBS, ZBS) 
B=HY12-(HY34-HY12)*TERM9+HINCYP(XBS, YBS, ZBS) 

GO TO 1375 


1150 CONTINUE 

TERM9 = ( XBS-X( IX) ) / TERK4 

C TERM9 IS -SCALE FACTOR FOR (EX,HY,HZ) (IH,J,K) 

C 

C 

ESTORE(L,NPT)=EX34+(EX34-EX12)*TERM9+EINCXP(XBS,YBS,ZBS) 
A=HZ34+( HZ 34 -HZ 12) *TER119+I1INCZP (XBS,YBS,ZBS) 
B=HY34+(HY34-HY12)*TERM9+HINCYP(XBS, YBS,ZBS) 

GO TO 1375 
C 
C 

C NPLANE = 3 OR 4 — XZ PLANES 

C 

C COMPUTE TERMS FOR CHARGE 

C 

C 

C DETERMINE IM,IX,KM,KX FOR EY 

C 

1155 CONTINUE 

DO 1160 1=1, NX1 

IF(XOBS(NPT).LT.XO(I+1)) GO TO 1165 
1160 CONTINUE 
1165 IM = I 

IX = 1+1 
DO 1170 K=1 ,NZ1 

IF(Z0BS(NPT).LT.Z0(K+1)) GO TO 1175 
1170 CONTINUE 
1175 KM = K 

KX = K+l 

IF(L.EQ.100)PRINT 410,IM,IX, JM, JX,KM,KX 
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XERMl = (ZBS - Z0(KM))/DZ0(KM) 

C TERMl IS SCALE FACTOR FOR EY(I,J,KX) 

TERM2 = (YBS - Y0( JM) ) /DY0( JM) 

TER113 = (XBS - X0(IM))/DX0(IM) 

C TERM3 IS SCALE FACTOR FOR EY(IX,J,K) 

TER1-14 = 0. 5*(DX0(IX)+DX0(IM)) 

TERM5 = 0. 5*(DY0( JX)+DY0( Jfl) ) 

TERM6 = 0. 5*(DZ0(KX)+DZ0(KM) ) 

EY1=EY(IM,JM,KM)+(EY(IM,JM,KX)-EY(IM,JM,KM))*TERM1 
EY2=EY( IX, JM,KM)+(EY( IX, JM,KX)-EY( IX, JM,KM) )*TERM1 
EY3=EY(IM,JX,KM)+(EY(IM,JX,KX)-EY(IM,JX,KM))*TERM1 
EY4=EY(IX,JX,KM)+(EY(IX, JX,KX)-EY(IX, JX,KM) )*TERM1 
C 

EY12 = EY1 + (EY2-EY1 )*TERM3 
EY34 = EY3 + (EY4-EY3)*TERM3 
C 

C DEBUG PRINT 

C 

IF(L.NE.IOO) GO TO 1202 

PRINT 500 , ( ( ( EY( I , J , K.) , I-IM , IX) , J-JM , JX) , K=KM,KX) 
PRINT 501, EY1 ,EY2,EY3 ,EY4 ,EY12,EY34 
1202 CONTINUE 

COMPUTE TERMS FOR CURRENT DENSITY 


DETERMINE KM,KX FOR HX IM,IX SAME AS FOR EY 
DO 1210 K=1 ,NZ1 

IF(ZOBS(NPT) .LT.Z(1I+1) ) GO TO 1220 
1210 CONTINUE 
1220 KM = K 
KX = K+l 

IF(L.EQ.IOO) PRINT 410,IM,IX,JM,JX,K1-1,KX 
TERMl = (ZBS - ZO(KM) ) /DZO(KM) 

TERM2 = (YBS - YO(JM) )/DY0( JM) 

TERM3 = (XBS - X0( IM) ) /DX0( IM) 

TERM3 IS SCALE FACTOR FOR HX(IX,J,K) 

TERM4 = 0. 5*(DX0( IX)+DX0( IM)) 

TERM5 = 0. 5*(DY0( JX)+DY0( JM) ) 

TERM6 = 0. 5*(DZ0(KX)+DZ0(KM) ) 

TERM7 = (ZBS-Z(KM) )/TERM6 

TERM7 IS SCALE FACTOR FOR HX(I,J,KX) 
liXl-HX ( IM , JM , KM) + ( IIX ( IM , Jt i , KX) - HX ( IM , JM , KM ) ) *TERM7 
11X2 =HX( IX , JM , KM)+( HX( IX , JM , KX) -11X( IX, JM , KM) ) *TER117 
HX3=HX( IM, JX,KM)+(11X(IM, JX,KX)-HX(IM, JX,KM) )*TERM7 
HX4=HX( IX, JX,KM)+(UX( IX, JX,KX)-HX( IX, JX,KM) )*TERM7 

HX12 = HX1 + (HX2-UX1)*TERM3 
11X34 = 11X3 + ( HX4-HX3 ) *TERM3 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1222 

PRINT 502 , ( ( ( HX( I , J , K) , 1= IM , IX) , J= JM , JX) , K=KM , KX ) 
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PRINT 503, HX1 , HX2 , HX3 , HX4 , HX1 2 , HX34 
1222 CONTINUE 


C DETERMINE IH, IX,KM,KX FOR HZ 

C 

DO 1225 1=1, NX1 

IF(XOBS(NPT).LT.X(I+l)) GO TO 1230 
1225 CONTINUE 
1230 IM = I 
IX = 1+1 
DO 1235 K=1 ,NZ1 

IF(|zOBS(NPT) .LT.Z0(K+1)) GO TO 1240 
1235 CONTINUE 
1240 KM = K 
KX = K+l 

IF(L.EQ.IOO) PRINT 410, IM, IX, JM, JX,KM,KX 
TERM1 = ( ZBS - ZO(KM))/DZO(KM) 

C TERM1 IS SCALE FACTOR FOR HZ(I,J,KX) 

TERM2 = ( YBS - YO(JM))/DYO( JM) 

TERM3 = (XBS - XO(IM))/DXO(IM) 

TERM4 = 0. 5*(DX0( IX)+DXO( IM) ) 

TERM5 = 0.5*(DY0( JX)+DYO(JM) ) 

TERM6 = 0.5*(DZ0(KX)+DZ0(KM) ) 

HZ 1=HZ ( IM , JM , KM)+( 11Z( IM , JM , KX)-HZ ( IM , JM , KM) ) *TERM1 
HZ2=HZ( IX, JM,KM)+(HZ( IX, JM,KX)-HZ( IX, JM,KM))*TERM1 
HZ3=HZ(IM,JX,KM)+(I1Z(IM,JX,KX)-11Z(IM,JX,KM))*TERM1 
HZ4=HZ( IX, JX,KM)+(HZ( IX, JX,KX)-HZ(IX, JX,KM))*TERM1 
C 

TERM8 = (XBS-X(IM) )/TERM4 
C TERM8 IS SCALE FACTOR FOR HZ(IX,J,K) 

HZ 12 = HZ1 + (11Z2-I1Z1)*TERM8 
HZ34 = 11Z3 + (HZ4-HZ3)*TERM8 
C 

C DEBUG PRINT 

C 

IF(L.NE.IOO) GO TO 1242 

PRINT 504,(((HZ(I,J,K) ,I=IM,IX) , J=JM, JX) ,K=KM,KX) 

PRINT 505, HZ1 ,HZ2 ,HZ3,HZ4,HZ12 ,HZ34 
1242 CONTINUE 
C 

IF(NPLANE(NPT).EQ.4) GO TO 1250 
C 

TERM9 = (Y(JM)-YBS)/TERM5 

C TERM9 IS -SCALE FACTOR FOR (EY,HZ,HX)(I,JX,K) 

C 

PRINT 995 

995 F0RMAT(*CHK1*) 

C 

C 

PRINT 996 

996 FORMAT (*CHK2*) 

C 

E STORE (L,NPT)=EY 12 - (EY34-EY12)*TERM9 + EINCYP(XBS,YBS,ZBS) 
PRINT 997 
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997 FORMAT ( *CHK3* ) 


A = HZ12 - (UZ34-HZ12)*TERM9 + HINCZP(XBS,YBS,ZBS) 
PRINT 998 

998 FORMAT (*CHK4*) 

B = HX12 - (11X34-11X12 )*TERM9 + 11INCXP(XBS,YBS,ZBS) 
PRINT 999 

999 FORMAT ( *CHK5* ) 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1375 
EINC = EINCYP(XBS,YBS,ZBS) 

111 INC = HINCXP(XBS,YBS,ZBS) 

112 INC = 11INCZP(XBS,YBS,ZBS) 

PRINT 513, EINC, H1INC,H2 INC 
PRINT 506,ESTORE(L,NPT) ,A,B 
GO TO 1375 


1250 TERM9 = (YBS-Y ( JX) )/TERM5 

TERM9 IS -SCALE FACTOR FOR (EY,HZ,HX) ( I, JM,K) 


EST0RE(L, NPT)=EY34 + (EY34-EY12)*TERM9 + EINCYP(XBS, YBS, ZBS) 
A = HZ34 + ( HZ34-HZ 12 ) *TERM9 + IIINCZP(XBS, YBS, ZBS) 

B = 11X34 + (11X34-11X1 2 )*TERM9 + liINCXP(XBS,YBS, ZBS) 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1375 
EINC = EINCYP(XBS , YBS, ZBS) 

1I1INC = HINCXP(XBS, YBS, ZBS) 

H2INC = HINCZP(XBS,YBS,ZBS) 

PRINT 513, EINC, HI INC, H2 INC 
PRINT 506 ,ESTORE(L,NPT) ,A,B 
GO TO 1375 
C 
C 

C NPLANE - 5 OR 6 — XY PLANES 

C 

C COMPUTE TERMS FOR CHARGE 

C 

C 

C DETERMINE IU,IX,JM,JX FOR EZ 

C 

1270 CONTINUE 

DO 1275 1=1, NX1 

IF(XOBS(NPT).LT.XO(I+1)) GO TO 1280 
1275 CONTINUE 
1280 IM = I 

IX = 1+1 
DO 1285 J=1 ,NY1 

IF( YOBS(NPT) .LT. Y0( J+l ) ) GO TO 1290 
1285 CONTINUE 
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1290 JM=J 

JX = J+l 

IF(L.EQ.IOO) PRINT 410,IM,IX, JM,JX,KM,KX 
TERM1 = ( ZBS - Z0(KM))/DZ0(KM) 

TERM2 - (YBS - Y0( JM) )/DY0( JM) 

C TERM2 IS SCALE FACTOR FOR EZ(I,JX,K) 

TERM3 = (XBS - X0(IM))/DX0(IM) 

C TERM3 IS SCALE FACTOR FOR EZ(IX,J,K) 

TERM4 = 0.5* (DX0( IX)+DX0( IM) ) 

TERM5 = 0. 5*(DY0( JX)+DY0( JM) ) 

TERM6 = 0.5*(DZ0(KX)+DZ0(KM) ) 

1300 EZ1=EZ( III, JM,KM)+( EZ( IX, JM, KM)-EZ( IM, JM ,KM) )*TERM3 
EZ2=EZ(IM,JX,KM)+(EZ(IX,JX,KM)-EZ(IM,JX,KM))*TERM3 
EZ3=EZ( IM, JM, KX)+( EZ( IX, JM,KX)-EZ( IM, JM, KX) )*TERI13 
EZ4=EZ( IM, JX,KX)+(EZ( IX, JX,KX)-EZ( IM, JX,KX) )*TERM3 
C 

EZ12 = EZ1 + (EZ2-EZ1 )*TERM2 
EZ34 = EZ3 + (EZ4-EZ3)*TERM2 


DEBUG PRINT 

IF(L.NE.IOO) GO TO 1302 

PRINT 509,(((EZ(I,J,K) ,I=IM,IX) , J=JM, JX) ,K=KM,KX) 
PRINT 510, EZ1 ,EZ2 ,EZ3 ,EZ4 ,EZ12,EZ34 
1302 CONTINUE 

COMPUTE TERMS FOR CURRENT DENSITY 


DETERMINE JM,JX FOR HX , IM,IX SAME AS FOR EZ 
DO 1310 J=1,NY1 

IF( YOBS(NPT) .LT .Y( J+l) ) GO TO 1315 
1310 CONTINUE 
1315 JM = J 
JX = J+l 

IF(L.EQ.IOO) PRINT 4 1 0 , IM , IX , JM , JX , KM , KX 
TERM1 = (ZBS - Z0(KM))/DZ0(KM) 

TERM2 = (YBS - Y0( JM) )/DY0(JM) 

TERM3 = (XBS - X0( IM) ) /DX0( IM) 

TERM3 IS SCALE FACTOR FOR HX(IX,J,K) 

TERM4 » U. 5* (DX0( IX)+DX0( IM) ) 

TERM5 = 0. 5* (DY0( JX)+DY0( JM) ) 

TERM6 = 0. 5*(DZ0(KX)+DZ0(KM) ) 

HX1=HX( IM , JM , KM)+( UX( IX , JM , KM)-ilX( IM , JM , KM) ) *TERM3 
HX2=HX(IM,JX,K11)+(HX(IX,JX,KM)-HX(IM,JX,KM))*TERM3 
1IX3=HX( IM , JM ,KX)+(11X( IX , JM , KX)-HX( IM , JM , KX) )*TERil3 
11X4=HX(IM,JX,KX)+(HX(IX,JX,KX)-UX(IM,JX,KX))*TERM3 

TERM7 = (YBS-Y(JM) )/TERM5 

TERM7 IS SCALE FACTOR FOR HX(I,JX,K) 

11X12 = HX1 + ( HX2-UX 1 ) *TERM7 
11X34 = 11X3 + ( 11X4-11X3 )*TER1-I7 

DEBUG PRINT 
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IF(L.NE.IOO) GO TO 1318 

PRINT 302 , ( ( ( UX( I , J , K) , 1=IM , IX) , J= JM , JX) , K=KM , KX) 
PRINT 503, HXi , HX2 , HX3 , 11X4 , HX1 2 , 11X34 
1318 CONTINUE 
C 
C 

C DETERMINE IH,IX,JM,JX FOR HY 

C 

DO 1320 1=1, NX1 

IF(XOBS(NPT) .LT.X( 1+1) ) GO TO 1325 
1320 CONTINUE 
1325 IM= I 

IX = 1+1 
DO 1330 J-l.HYl 

IF(YOBS(NPT) . LT.Y0( J+l) ) GO TO 1335 
1330 CONTINUE 
1335 JM = J 
JX = J+l 

IF(L.EQ.IOO) PRINT 410, 1M,IX, JM, JX,KM,KX 
TERM1 = (ZBS - Z0(KM))/DZ0(KM) 

TERM2 = (YBS - Y0( JM) )/DY0(JM) 

C TERM2 IS SCALE FACTOR FOR HY(I,JX,K) 

TERM3 = (XBS - X0( IM) ) /DXO( IM) 

TERM4 = 0.5*(DX0(IX)+DX0(Iil) ) 

TERM5 = 0 . 5*(DY0( JX)+DY0( JM) ) 

TERM6 - 0. 5*(DZ0(KX)+DZ0(KM) ) 

TERM8 = (XBS— X( IM) ) /TERM4 
C TERMS IS SCALE FACTOR FOR HY(IX,J,K) 

HY1=HY( IM, JM,KM)+(UY( IX, JM,KM)— 11Y( IM, JM.KM) )*TERM8 
HY2=liY ( IM , JX , KM)+( HY ( IX , JX , KM) — 11Y ( IM , JX , KM) ) *TERM8 
11Y3=HY( IM , JM , KX)+( UY( IX, JM , KX)-HY( IM , JM, KX) )*TERM8 
HY4=HY( IM, JX,KX)+(HY( IX, JX,KX)— HY( IM, JX,KX) )*TERM8 
C 

11Y1 2 = HY1 + (11Y2-11Y1 )*TERM2 
hY34 = HY3 + (HY4-1IY3)*TERM2 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1340 

PRINT 511,(((UY(I,J,K) ,I=IM,IX) ,J=JM,JX) ,K=KM,KX) 
PRINT 512, HY1,HY2,HY3,HY4,HY12,HY34 
1340 CONTINUE 

IF(NPLANE(NPT) .EQ.6) GO TO 1350 

TERM9 = (Z(KM)-ZBS)/TERM6 

TERM9 IS -SCALE FACTOR FOR (EZ,HX,HY)(I, J,KX) 


EST0RE(L, NPT)=EZ12 - (EZ34-EZ12)*TERM9 + EINCZP(XBS,YBS,ZBS) 
A = 11X12 - (ilX34-HX12)*TERM9 + 1IINCXP(XBS, YBS, ZBS) 

B = HY12 - (11Y34-HY12 )*TERM9 + 11INCYP(XBS,YBS,ZBS) 

DEBUG PRINT 
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c 

IF(L.NE.LOO) GO TO 1375 
EINC = EINCZP(XBS,YBS,ZBS) 
1I1INC= HINCXP(XBS, YBS.ZBS) 
H2INC= HINCYP(XBS, YBS.ZBS) 
PRINT 513 , EINC ,111 INC ,H2 INC 
PRINT 506 ,ESTORE(L,NPT) ,A,B 
GO TO 1375 


1350 TERM9 = (ZBS-Z(KX) ) /TERM6 

TERM9 IS -SCALE FACTOR FOR (EZ , HX ,HY) ( I , J , KM) 


C 


C 


ESTORE(L,NPT)=EZ34 + (EZ34-EZ12)*TERM9 + EINCZP(XBS, YBS, ZBS) 
A = 11X34 + ( 11X34-11X1 2 ) *TERM9 + HINCXP(XBS, YBS.ZBS) 

B = t!Y34 + (UY34-HY12)*TERM9 + HINCYP(XBS,YBS,ZBS) 

DEBUG PRINT 

IF(L.NE.IOO) GO TO 1375 
EINC = EINCZP( XBS , YBS , ZBS ) 

U1INC= HINCXP( XBS, YBS.ZBS) 

112INC= HINCYP( XBS , YBS , ZBS) 

PRINT 513, EINC, H1INC.H2INC 
PRINT 506 ,ESTORE( L, NPT) ,A,B 
GO TO 1375 


1375 CONTINUE 

SINAN = SIN (ANGLE) 

COSAN = COS( ANGLE) 

HSTORl(L.NPT) » A*COSAN+B*SINAN 
HST0R2(L,NPT) = B*COSAN-A*SINAN 
CUMM( NPT )=CUMM( NPT )+HSTORl ( L, NPT) 
C 


IF(L.EQ.l) GO TO 1390 
IF(L.NE.IOO) GO TO 1400 
1390 CONTINUE 

PRINT 507, 11ST0R1(L,NPT),HST0R2(L,NPT) ,ESTORE(L,NPT) .ANGLE 
PRINT 508, SINAN, COSAN 
1400 CONTINUE 


C 

C FORMAT STATEMENTS 

C 

401 FORMAT ( 15) 

402 FORMAT (215, 4F 10. 3) 

407 F0RMAT(7H1 XOBS,/ ,3(10(F8.2,4X) ,/)) 

408 FORMAT (7 UO YOBS ,/ , 3( 10( F8 . 2 ,4X) , /) ) 

409 FORMAT (7 HO ZOBS,/ ,3(10(F8.2,4X) ,/)) 

410 FORMAT (1 OHO IM IX,8X,7liJM JX,8X,7HKM KX, /, 3(215 ,5X) ) 

500 FORMAT (6110 EY ,8E12.4,/) 

501 FORMAT ( 1110, 4X, 3HEY1 , 9X, 311EY2 ,9X, 3HEY3 ,9X, 3HEY4.8X, 4HEY12 ,8X, 

1 4HEY34 , / , 6E 1 2 . 4 , / ) 

502 FORMAT (6110 HX ,8E12.4,/) 

503 FORMAT ( 1110 , 4X, 3HHX1 , 9X, 31111X2 , 9X.3HHX3 , 9X, 3IHIX4, 8X, 4HHX12 , 8X, 
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1 41111X34 , / , 6E12 .4 , / ) 

504 FORMAT ( OHO HZ .8E12.4,/) 

505 FORMAT ( 1H0 , 4X ,31111X1 , 9X , 31U1Z2 , 9X , 3H11Z3 , 9X , 311H24 , 8X, 4HHZ12 , 8X, 

1 411HZ34,/ , 6E12 .4 , / ) 

506 FORMAT (9 HO E,H1 ,H2 , 5X, 3E12 . 4) 

507 FORMAT ( 16H0 111 ,H2 ,E, ANGLE, , 5X, 4E12 .4) 

508 F0RI'1AT( 14110 SIN (ANGLE) =,F12.3,5X, 12HC0S(ANGLE) =,F12.3) 

509 F0RMAT(6110 EZ ,8E12.4,/) 

510 FORMAT ( 1H0 ,4X, 3HEZ1 , 9X, 3HEZ2 ,9X,3HEZ3 , 9X, 31IEZ4,8X,4I1EZ12 ,8X, 

1 4HEZ34,/ ,6E12.4,/) 

511 FORMAT (6 HO 1IY ,8E12.4,/) 

512 FORMAT ( 1H0,4X, 3HHY1 , 9X, 3ilHY2 , 9X, 3HHY3.9X, 3HHY4,8X, 4H1IY12,8X, 

1 41IIIY34,/ ,6E12.4,/) 

513 FORMAT ( 1110 , 8H EINGP =,E12.4,2X,8HH1INCP =,E12.4,2X, 

1 8I1H2INCP =,E12.4) 

514 FORMAT ( 6H0 EX ,8E12.4,/) 

515 FORMAT ( 1H0 , 4X , 3HEY1 , 9X , 31IEY2 , 9X , 31IEY3 , 9X , 3HEY4 , 8X, 4HEY 12 , 8X , 

1 4HEY34 , / , 6E12 . 4 ,/) 

RETURN 

END 
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3.5.11 Subroutine PRINOUT 


PRINOUT prints out the data saved in DATASAV. It also extends the 

data record for H-field loops that define wire currents. The data is summed 

over time, which averages out any resonances (these behave as A^e “n^sin 1 ^ and 

when summed over several cycles yield approximately zero.) This leaves the sum 

of the decaying exponential response A e ^^dt, which is S = / T A e °o t dt. The 

o oo 

decaying exponential response is the portion of the total response that tracks 

the driving function which decays as e °o t . Thus, a Q is known from the 

excitation source and A can be found from S, the numerically derived sum, and 

T p, ^ 

the expression S = / A e o dt. With A and a known, the data record is then 

o o oo 

extended in time. This allows transforms to be made for data records as long 
as 50 ysec or a lower frequency limit of 10 KHz in the spectrum. 
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SUBROUTINE PRINOUT 


C 

COMMON/ EXTRAS/NX, NY ,NZ ,NX1 ,NY1 ,NZ1 ,N , M,MQ ,DT, XMU ,EPSO, EPS , NPTLIM , 
1 NN,NPTS, LMAX, SIGMA, C,T, PI, EXPF AC, IP, TX,TY,TZ, AMP, ALPHA, BETA, IDLS 
COMMON/ GRID/ X( 2 8) ,Y(28) ,Z(28) ,X0(29) ,YO(29) ,ZO(29) , 

1 DX(29) ,DY(29) ,DZ(29) ,DXO(28) ,DY0(28) ,DZ0(28) , 

2 DXI(29) ,DYI(29) ,DZI(29) ,DX0I(28) ,DY0I(28) ,DZOI(28) 

COMMON / OUT/ ESTORE (250,24) ,HSTOR1(2500 ,24) ,HSTOR2(250,24) , 

1 TST0RE(250) 

COMMON/ CUMMUL/ CUMM(24) ,LOC(24) 

COMMON/ PERM/ MM 
DIMENSION CSTORE(2500,8) 

REAL PP 
C 
C 

C PRINT OUTPUT 

C 

DO 1500 NPT = 1 ,NPTS 
C 

L=N/IP 

PP=IP 

C 

PRINT 24 ,LOC( NPT) 

24 FORMAT (8H1 LOC = ,13) 

PRINT 27 ,L 

27 FORMAT (80X/36X, 14 ,40X) 

PRINT 23, (TSTORE(LL) ,LL=1 ,L,3) 

23 FORMAT( IX, 1P6E13. 6 , IX) 

PRINT 27, L 

PRINT 23, (ESTORE(LL.NPT) ,LL=1 ,L,3) 

PRINT 27, L 

PRINT 23, ( HSTOR1 (LL.NPT) ,LL=1 ,L , 3 ) 

PRINT 27, L 

PRINT 23, (HSTOR2(LL, NPT ) ,LL=1 ,L, 3) 

C 

TT=T 

C 

DO 1501 LL=L, LMAX 
TT=TT+DT*8. 

HSTORI (LL , NPT ) =CUMM( NPT ) *DT*8 . * 

1 ( 1 . 0/ ( ( 1 .0-EXP ( -ALPdA*T) ) /ALPHA- ( 1 . 0-EXP ( -BETA*T) ) / BETA) ) * 

2 EXP ( -ALPHA*TT ) 

1501 CONTINUE 

C 

1500 CONTINUE 
C 

IF(M.NE.l.OR.MM.NE.l) GO TO 1600 
C 

DO 1510 LL=1 , LMAX 
CSTORE ( LL , 1 ) =HST0R1 ( LL , 1 ) 

C STORE (LL, 2)=HST0Rl (LL, 2) 

1510 CONTINUE 
C 

PRINT 24 ,LOC( 1) 

PRINT 23, (CSTORE ( LL , 1 ) , LL= 1 , LMAX , 3 ) 
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PRINT 24,LOC(2) 

PRINT 2 3 , ( CSTORE ( LL , 2 ) ,LL=1 ,LMAX,3) 
C 

GO TO 500 
C 

1600 CONTINUE 
C 

IF(M.NE.2.0R.MM.NE. 12) GO TO 1700 
C 

DO 1610 LL=1 , LIT AX 

CSTORE ( LL , 1 ) =HST0R1 ( LL , 1 ) 

CST0RE(LL,2)=HST0R1(LL,2) 

1610 CONTINUE 
C 

PRINT 24 ,LOC( 1 ) 

PRINT 23 , (CSTORE ( LL, 1) ,LL=1 ,LMAX, 3 ) 

PRINT 24 ,L0C(2) 

PRINT 23 , ( CSTORE ( LL ,2) ,LL=1 ,LMAX,3) 
C 

GO TO 500 
G 

1700 CONTINUE 
C 

I F ( M . NE . 2 ) GO TO 1800 
C 

DO 1710 LL= 1 , LilAX 

CST0RE(LL,1)=HST0R1(LL,1) 

CSTORE ( LL , 2 )=HST0R1 ( LL , 2 ) 

1710 CONTINUE 
C 

PRINT 24,LOC(l) 

PRINT 23 , ( CSTORE (LL, 1 ) ,LL=1 , L11AX, 3 ) 

PRINT 24,LOC(2) 

PRINT 23 , (CSTORE (LL, 2) ,LL=1 ,LMAX,3) 
C 

GO TO 500 
C 

1800 CONTINUE 


500 CONTINUE 
C 

RETURN 

END 

/EOR 
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APPENDIX A 


The governing equations for the lossy dielectric finite difference 
code are in the form 


e 



+ 



/ 3H 
\ 9y 


° 9H 

z l 

9z 



(e-e o ) 



where the x component of the E-field is used, for example. Consider now a 
plate of fractional thickness T in the x direction, i.e. the plate is TAx 
thick, so that T is a fraction of the cell dimension in the desired direction. 
The effective difference equation for the cell in which the plate is present 
can be found from the average of the above difference equation for a region in 
which the plate is present (a=0, e=e o ) and where there is only free space (o= 0, 
e=eQ> . More precisely 



In explicit finite difference form this is expressed as 
E®(I,J,K) n+1/2 = E^(I,J,K) n “ 1/2 e" oTAt/£eff 

+ O-e" 0TAt/£ eff )x 

|"-E^(I,J,K) n ~ E^(I,J,K) n 

H®(I,J,K) n - H®(I,J-1 ,K) n 
+ oT( Y( J) - Y(J-1)) 


Hy(I , J, K) n - H®( 

— JLi J 


oT(Z(K) - Z( 


I , J , K-1 ) n “ 
K-D) 


where e eff = C ( 1 — T ) e Q + Te]. 


Similar expressions hold in the outer directions for the thin plate formalism. 


139 


REFERENCES 


1. Yee, K. S. , "Numerical Solution of Initial Boundary Value Problems 
Involving Maxwell's Equations in Isotropic Media," IEEE Trans. Ant. & 
Prop. , Vol . AP-1 4 , pp. 302-307, May 1966. 

2. Merewether, D. E. , "Transient Currents Induced on a Metallic Body of 
Revolution by an Electromagnetic Pulse," IEEE Trans. EMC , Vol. EMC-1 3 » PP* 
41-44, May 1971 . 

3. Taylor, C. D. ; Lam, D. H. and Shumpert, T. H. , "Electromagnetic Pulse 
Scattering in Time-Varying Inhomogeneous Media," IEEE Trans. Ant. Prop. , 
Vol. AP-1 7 , No. 5, September 1969. 

4. Lax, P. D., "Differential Equations, Difference Equations and Matrix 
Theory," Commun. Pure Appl. Math. , Vol. 2. pp. 175-194* November 1958. 

5. Holland, R., "THREDE: A Free-Field EMP Coupling and Scattering Code," 

Mission Research Corporation, AMRC-R- 85 , September 1976. 

6. Kunz , K. S. and Lee, K. M. , A Three-Dimensional Finite-Difference Solution 

of the External Response of an Aircraft to a Complex Transient EM 
Environment: Part 1 - The Method and Its Implementation; and Part 2 - 

Comparison of Predictions and Measurements," IEEE Transactions EMC , Vol. 
EMC-20, No. 2, May 1978. 


140 



1. Report Na 


2 . 


Government Accession No. 


3. 


Recipient's Catalog No. 


4. 


NASA CR-166079 | 

Title and Subtitle 

Generalized Three-Dimensional Experimental Lightning 
Code (G3DXL) User’s Manual 


5. Report Date 

February 1986 

6. Performing Organisation Code 


7. Author(s) 


8 . 


Performing Organisation Report No. 


Karl S . Kunz 


9. Per forming Organisation Name and Address 

Kunz Associates Incorporated 
P.0. Box 14373 
Albuquerque, NM 87191 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 

15. Supplementary Notes 


KAI-R-1 

10. Work Unit No. 


11. Contract or Grant No. 

NAS1-16591 

13. Type of Report and Period Covered 

Contractor Report 

1981-1986 

14. Sponsoring Agency Code 

505-66-21-04 


Langley Technical Monitor: Felix L. Pitts 


16. Abstract 

A new electromagnetic coupling computer code is an integration of several 
finite-difference codes based on THREDE, an established electromagnetic 
computer modeling technique. For a given lightning model, the new integrated 
code will be capable of predicting, with better resolution than was hitherto 
available, EM fields caused by lightning at several locations outside or 
inside metal or composite aircraft. Increased resolution is obtained by 
making two computer runs: total (unexpanded) geometry and expanded 

geometry. The total response is found and then is expanded, allowing 
greater resolution in geometry and higher frequency response. 


17. Key Words (Suggested by Author(s)) 


18 . 


Distribution Statement 


Finite difference code. Lightning, 
LEMP, Computer modeling 



Subject category - 47 


19. Security Qesif. (of this report) 

Unclassified 


20. Security Oae*f . (of thie peg*) 

Unclassified 


21. No. of Pages 

144 


22. P*»ce 


















